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an  assumed  "analytical"  transmittance  function,  using  the 
same  type  of  data.  Computerized  numerical  techniques  are 
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a  secondary  effort  a  structural  breakdown  of  the  Lowtran 
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I.  Introduction 

Following  the  efforts  of  Elssaser^  numerous  workers 
have  attempted  to  arrive  at  computationally-simple  models 
for  gaseous  molecular  transmittance,  averaged  over  narrow 
spectral  '  intervals  in  the  infrared.  These  efforts  may  be 
naturally  divided  into  those  involving  the  analytical 
derivation  of  a  mean  transmittance  function  from  Beer's 
Law,  and  those  involving  the  extraction  of  the  transmittance 
function  itself  from  transmittance  data.  Traditionally, 
the  former  are  called  "analytical"  and  the  latter  are  called 
"empirical".  The  method  normally  used  in  the  empirical 
models  consists  of  the  extraction  of  the  transmittance 
function  through  graphical  techniques,  with  the  adoption 
of  a  relationship  between  spectral  and  absorber  parameters. 
In  the  development  of  analytical  models  a  transmittance 
function  is  adopted  at  the  offset,  and  the  spectral  and 
absorber  parameters  are  afterward  determined  through  com¬ 
puterized  numerical  procedures. 

In  the  work  reported  here  the  authors  present  a 
totally  computerized  version  of  the  classical  graphical 
methods  for  the  extraction  of  the  empirical  transmittance 
function.  This  is  followed  by  a  presentation  of  a  numerical 
method  which  uses  a  double-exponential  transmittance  func¬ 
tion  for  the  development  of  analytical  band  models.  Both 
methods  are  then  applied  to  20  cm  ^  averaged  line-by-line 


2 


transmittance  data  for  the  atmospheric  trace  gases  SC^, 

NO,  NC^,  and  NH^ .  The  model  parameters  are  listed  at 
5  cm  1  intervals  throughout  the  major  absorption  bands  of 
these  gases  for  the  convenience  of  the  community  of  band 
model  users.  Although  the  methodology  is  applied  speci¬ 
fically  to  the  trace  gases,  no  restrictions  are  immediately 
evident  in  the  extension  to  other  gaseous  absorbers  in  the 

infrared.  In  fact,  the  analytical  method  was  successfully 
2 

applied  earlier  to  the  principal  band  centers  of  the  major 
absorbers  vapor,  0^  and  the  uniformly-mixed  gases. 

As  an  application  of  the  results  found  through  this 
effort,  the  band  models  for  the  trace  gases  were  incor¬ 
porated  in  the  widely-used  code  called  Lowtran.  To 
facilitate  the  inclusion  of  these  models,  as  well  as  of 
others,  the  code  was  broken  down  into  separate  subroutines 
or  modules  controlled  by  a  master  program.  The  subroutines 
include  the  evaluation  of  the  equivalent  absorber  amount, 
the  selection  of  the  sp ec t r a  1 ly-e f f e c t ive  attenuation  model 
and  the  individual  attenuation  models.  The  principal 
purpose  of  the  modularization  is  to  assist  users  with  the 
modification  of  the  code  to  suit  their  individual  re¬ 
quirements  on  transmission  models. 
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II.  The  Transmittance  Equation 

The  monochromatic  transmittance  at  frequency  V 
for  the  passage  of  infrared  radiation  through  a  path  length 
Z  in  an  inhomogeneous  medium  with  pressure  and  temperature 
distributions  P(Z)  and  T(Z),  respectively,  is  given  by 
Beer's  Law  in  the  form 

Tu  -  .  -j  VP-T>dt,<Z>  (1) 

where  is  the  resultant  absorption  coefficient  for  all 
contributing  lines  and  gaseous  absorbers,  and  U  is  the 
absorber  amount.  For  broadband  radiation  detected  by  an 
instrument  of  spectral  response  <j>^ ,  the  variable  of  interest 
is  the  weighted  mean  transmittance  t,  defined  as 


T 


T  <f>  dv  / 

vvv 


<j>vdv 


(2) 


Equation  (2)  has  been  evaluated  analytically  over  a 

spectral  interval  Av  for  the  special  case  of  Lorentzian 

broadened  lines  having  assumed  line  distributions  and 

1  3 

intensities,  leading  to  the  classical  band  models  ’  . 

Numerous  variations  of  the  classical  band  models  may  be 

found  in  the  literature,  most  of  which  specify  the  analytical 

form  of  x  in  terms  of  mean  line  or  meteorological  variables. 

4 

A  notable  exception  is  the  model  of  King  which  expresses 
the  homogeneous-path  transmittance  as 

T  =  g ( SanU ) , 


(3) 


where  g  is  a  function  to  be  determined  empirically,  S  is 
the  mean  line  intensity,  a  is  the  mean  line  half-width  and 
n  is  an  absorber  parameter  with  the  physical  constraints 
of  zero  and  one  in  the  weak-line  and  strong-line  limits, 
respectively.  The  path  inhomogeneity  may  be  accounted  for 
in  Eq.  (3)  through  the  Curtis-Godson  equivalences 

SanU  =  |  S(Z)ctn(Z)dU(Z)  .  (4) 

From  practical  considerations,  it  is  often  desirable 
to  transform  the  argument  in  Eq.(3)  with  the  known  relations 


(^) 

(5) 

T  h 

(~>  (y2) 

(6) 

o 


in  order  to  obtain 


n  T  m 

T  =  g  {  C(~ )  (^)  U  }  ,  (7) 

o 

where  C  is  a  spectral  parameter  combining  Sq  and  a”,  m  is 
an  absorber  parameter  combining  the  temperature  exponents 
of  S  and  a,  a  is  an  absorber  constant,  and  the  subscript 
"o"  denotes  standard  conditions.  For  computational  con¬ 


venience  Eq.  (7)  may  be  expressed  as 


where 


x  -  Cf  +  log1QW 
C-  >  log10C 


(9) 

(10) 

(ID 


Here,  f  is  the  transmittance  function,  C'  is  the  spectral 
parameter,  W  is  the  equivalent  absorber  amount,  and  n  and 
m  are  the  absorber  parameters;  all  of  which  are  to  be 


determined  from  transmittance  data  for  each  absorber. 
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III.  Computerized  Method  of  Empirical  Model  Development 
3.1  Introduction 

Assuming  the  availability  of  equal  transmittance 
data,  which  is  defined  below,  we  have  developed  an  algorithm, 
called  ADSET,  which  evaluates  absorber  parameters  n,  m, 
spectral  parameters  C'(v)  and  an  empirical  transmission 
function  simultaneously.  In  the  algorithm  the  transmission 
function  is  linearized  and  a  linear  regression  technique 
is  utilized  for  parameter  evaluation.  In  order  to  evaluate 
the  band  model  parameters  and  the  empirical  transmission 
function  simultaneously,  a  set  of  auxiliary  variables  are 
introduced.  Each  data  point  is  identified  through  the 
auxiliary  variables  to  an  absorption  band  and  to  a  trans¬ 
mittance  ’cut1.  This  enables  us  to  obtain  globally  optimal 
set  of  parameters  and  the  empirical  transmission  function 
simultaneously. 

Based  on  the  derived  optimal  pointwise  transmission 
function,  a  piecewise  analytical  transmission  function  is 
developed.  The  commonly  used  computer  code  Lowtran  for 
the  evaluation  of  atmospheric  transmittance  can  be  greatly 
simplified  by  the  use  of  this  piecewise  analytical  trans¬ 
mission  function  to  model  the  major  absorbers. 

Finally,  the  code  ADSET  also  contains  a  subroutine 
which  can  compute  the  spectral  parameter  value  C'(v)  for 
non-major  absorption  bands. 


3.2 


Data  Structure 


Several  transmittance  values  x^  ,  j-1,  2,  ....  NCUT 
are  chosen  a  priori,  where  NCUT  is  the  number  of  chosen 
transmittance  values.  Curves  of  growth  data  (i.e.  T  versus 
U)  for  each  layer  of  atmosphere  are  assumed  to  be  given  at 
these  transmittance  values.  Therefore,  the  curves  of 
growth  have  'cut'  structure,  namely,  all  data  points  are 
on  one  of  the  transmittance  cuts  r  =  Xj  ,  j«l,  2,  ...,  NCUT 

(See  Fig.  1).  We  call  a  data  set  with  this  cut  structure 
an  'equal  transmittance'  data  in  the  sequel. 

3.3  Linearization  of  Transmission  Function 

Since  f  in  Eq.(8)  is  known  to  be  strictly  monotone 
decreasing  from  one  to  zero  as  x  changes  from  to  °°, 
there  exists  an  inverse  function  f  1  defined  on  (0,1)  such 
that 

x  =  f  _  1  (  x  ) 

=  C'  +  logW 

T 

=  C1  +  nlog(|“)  +  mlog(^)  +  logU.  (12) 

o 

Let  us  define  x^ ,  j=l,  2,  ...»  NCUT  be  the  inverse  image 

of  the  prechosen  transmittance  values  x^  ,  j’l,  2,  ..., 

NCUT  i.e., 

Xj  =  f"1(xj),  j=  1,  2 . NCUT,  (13) 


(LOGARITHMIC  SCALE) 


total  square  error  E  are  given  by 
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E  =  z 

j=l 


*ij 


Lii 


j=l  k=l 


'ijk* 


(17) 


NB  NB  J  i  Lij 

E  =  E  Ei  =  E  ,  ?  .  ?  .  Eijk‘ 


i=  1 


1=1  j=l  k=l 


(18) 


where  and  NB  are  the  numbers  of  the  cuts  in  i-th 

absorption  band  and  of  the  absorption  bands,  respectively, 

The  final  expression  can  be  simplified  if  we  assume  that 

the  number  of  layers  (L^  )  in  every  cut  is  equal  to  a 

constant  L .  .  For  this  case 
i 


NB  J  Li 


E  =  E 


i-1  j-1  k= 1 


E  .  ... 

ijk 


(19) 


Our  objective  is  to  find  optimum  set  of  parameters  (n*, 

XNCUT*>  WhiGh 

minimizes  this  grand  total  error  E. 


m*,  C^*,  C'2* 


r  '  * 

’  ^  NB  ! 


Xx*,  x2* 


3.5 


Auxiliary  Variables 


In  order  to  perform  the  minimization  of  E  with 

respect  to  the  above  parameters  simultaneously,  \.e  modify 

the  square  error  E^^  30  that  it  contains  all  the  parameters, 

This  is  done  by  introducing  two  sets  of  auxiliary  variables 

u^,  i*»l ,  2,  ...,  NB  and  v^  ,  j  =  l,  2,  NCUT.  Using  them, 

E...  is  redefined  as 
ijk 


'ijk 


=  {nlog(-j^)  +  mlogC^M  +  u1>ijk  C|  +  ...  +  uHB(ijkC* 

2 

+  Vl,ijkKl  +  +  VNCUT,ijk  KNCUT  “  ("l08Uijk)}  (20) 


NB 
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where  Kj  -  “Xj '  3“^*  2,  ....  NCUT,  The  auxiliary  variables  act  as 
identifiers  of  the  band  and  the  cut.  Tf  a  data  point  is  for  1-th 
band's  J-tb  cut,  then  u^  -  1  and  u^  =  0  for  all  1  4  1  and  v^  =  1 
and  Vj  -  0  for  all  j  4  5.  Thus,  only  the  spectral  parameter  C'^ 
and  the  cut  parameter  Kj  corresponding  to  the  current  data  are  active 
and  all  other  spectral  and  cut  parameters  disappear.  Hence,  Eq.  (20) 
reduces  to  Eq.  (15).  The  change  from  x.  to  K  =  ~'.<A  is  made  in  order 

J  J  J 

to  symmetrize  the  coefficient  matrix  of  the  resulting  normal  equation. 
This  change  makes  it  possible  to  utilize  any  specialized  solution 
method  for  the  symmetric  normal  equation  when  the  space  conservation 
is  important. 


3.6  Regression  Analysis 

Using  the  grand  total  error  E  with  the  redefined  in 

Eq.  (20),  the  best  parameter  values  n*,  m*,  C'-^*,  c'nb*»  k*»***» 

jlf 

KNCUT  are  slmultaneously  determined  by  the  linear  regression.  Setting 
the  partial  derivatives  of  E  with  respect  to  parameters  equal  to 
zero  results  in  a  linear  normal  equation  of  the  form  AX  "  B,  where 
A,  B  and  X  are,  respectively ,  a  symmetric  coefficient  matrix,  a 
constant  vector  and  a  parameter  vector  defined  by 


*  2 

•Eu  „ 


a  (log  ij£) 


"(log  f-)2 

o 
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B  =  jl(-vNCUT  logU),  EC-v^ogU),  I(-uNBlogU) . 

T  t 

Z(-u2logU),  K-log(^)logU),  Z(-log(|-)loguj  ,  (22) 

o 

X  =  ^NCUT*  *  *  *  *  ^1*  ^NB*  ^  2  ’  m  *  (23) 

The  *  in  Eq  .  (21)  represents  some  nonzero  elements.  Also 

the  Z  in  the  above  equations  represents  the  triple  sum 
NB  J  L. 

Z  Z  Z  in  Eq.  (19).  One  may  realize  that  C'  does  not 
i=l  j-1  k=l 

appear  in  Eq .  (23)  and  hence  the  corresponding  auxiliary 

variable  u^  is  also  absent  from  Eqs.  (21)  and  (22).  This 

is  because  one  of  Cl,  ...»  C '  is  dependent  on  other  C! 

1  N  B  i 

so  that  C' ,  ....  C'  cannot  be  determined  uniquely.  It  is 

1  N  B 

necessary  that  one  of  C^s  be  given  a  number  a  priori.  Here 
Cj  is  chosen  and  is  given  the  value  zero,  and  therefore  is 
eliminated  from  the  parameter  vector  X.  This  choice  calls 
for  some  explanation.  On  t  vs,  logW  diagram  the  optimum 
empirical  transmission  function  can  be  placed  anywhere. 

What  it  amounts  to  is  that  a  different  placement  results 
in  a  different  set  of  C|  values  which  is  a  linear  shift 
(addition  or  subtraction  of  a  constant)  of  another  set  of 
values.  Only  the  relative  relationship  among  is 
unique.  This  is  clearly  indicated  in  Fig.  2. 

Since  the  placement  of  the  empirical  transmission 
function  is  arbitrary,  we  may  position  it  on  the  data 
points  corresponding  to  the  first  absorption  band.  In 
other  words,  the  first  absorption  band  is  taken  as  the 


j 
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reference  band  and  the  corresponding  spectral  parameter 
Cj  is  set  to  be  zero. 

The  queuing  of  parameters  in  X  vector  is  determined  in 
such  a  way  that  as  many  upper  principal  minor  matrices  as 
possible  become  diagonal  (See  Eq .  (21)).  This  arrangement 

can  reduce  the  amount  of  computation  in  the  early  stage  of 
Gauss  elimination  steps  when  the  normal  equation  is  solved, 
and  can  result  in  less  computational  error. 

3.7  Piecewise  Analytical  Transmission  Function 

After  the  best  parameter  values  are  computed,  the 
piecewise  analytical  transmission  function  is  generated 
by  the  piecewise  interpolation.  The  transmittance  region 
(0,1)  is  divided  into  NCUT  -  1  subregions  by  the  trans¬ 
mittance  cuts  T  2  ,  T3,  ...,  Let  T1  >  T2  >  > 

TNCUT’  t^en  t'ie  subregions  are  given  by 

Subregion  1  [ T2 , 1  ) , 

Subregion  2  ^T3’X1^’ 

Subregion  NCUT-2  ^ TNCUT-1 ’ TNCUT-2 ^ ’ 

Subregion  NCUT-1  ’ XNCUT-1 ^ ' 

The  top  and  bottom  subregions  contain  t  ^  and 


as  an  inner  point,  respectively.  The  interpolation  in  each 
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subregion  is  done  by  the  double  exponential  function  defined 
by 

.  .  2 

a  +  a  ~x  +  a  .  x 

x(x)  =  exp  {-10  }  .  (24) 

The  generally-used  linear  interpolation  is  not  used  here 

since  subregions  cannot  be  assumed  small  enough  for  the 

linear  approximation  to  be  valid.  Furthermore,  the  linear 

interpolation  is  totally  inadequate  for  the  top  and  bottom 

subregions.  On  the  other  hand,  the  double  exponential 

function  takes  the  values  between  and  is  asymptotic  to  one 

and  zero  as  the  argument  varies  from  -°°  to  It  is  also 

known  that  this  function  closely  approximates  the  standard 

2  1 

empirical  transmission  function  used  in  the  Lowtran  code.  ’ 
The  parameters  a^ ,  a 2,  and  a^  for  each  subregion 
are  determined  by  two  different  methods.  The  first  method 
assumes  that  a^  =  0  and  uses  no  further  data  to  compute 
a^  and  .  They  are  simply  determined  by  the  condition  that 
the  interpolation  function  in  each  subregion  passes  through 
the  end  points.  In  the  top  and  bottom  subregions,  the  func¬ 
tion  is  required  to  pass  through  two  points;  (t^,  x^)  and 

(t2,  x  2  >  for  the  top  and  ^TNCUT-1,XNCUT-1'  and  ^NCUT’ 

XNCUT^  ^°r  ttie  bottom  subregions. 

The  second  method  does  not  assume  that  a^  -  0  and 
requires  additional  data  to  compute  parameter  values.  The 
same  condition  that  each  interpolation  function  passes 


17 


IV.  Computerized  Method  of  Analytical  Model  Development 

4.1  Introduction 

In  the  last  chapter,  we  assumed  no  analytical  form 
for  the  transmission  function  t  =  f  (x)  when  the  standard 
transmission  function  was  computed.  Here,  by  assuming  the 
double  exponential  form  given  by  Eq .  (24)  as  the  transmission 

function  for  the  entire  transmittance  range,  we  derive  an 
algorithm  which  can  evaluate  the  best  function  parameter 
values  a^,  and  a^;  together  with  the  band  model  para¬ 

meters  n,  m,  and  C|.  Note  that  the  double  exponential 
function  was  used  for  the  piecewise  interpolation  in  the 
last  chapter.  But  the  computation  of  the  function  para¬ 
meters  was  performed  after  the  band  model  parameters  and 
the  empirical  transmission  function  were  obtained.  In 
other  words,  the  computation  in  the  last  chapter  was 
sequential  but  not  simultaneous.  The  algorithm  we  present 
in  this  chapter  is,  on  the  contrary,  the  simultaneous 
evaluation  of  all  parameters.  The  preliminary  development 
of  this  algorithm  can  be  found  in  Ref.  5. 

4.2  Basic  Equations 

The  basic  equations  are  Eq .  (8)  and  Eq.  (24)  of 

the  last  chapter,  which  are  cited  here  for  the  ease  of 
reference . 

T  =  f(x)  2  (27) 

a  +  a„x  +  a.x 

f(x)=exp(-10  }. 


(28) 
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Now,  since  we  have  assumed  the  function  form,  we 
can  compute  the  transmittance  if  we  have  the  value  of  x. 
Hence,  we  do  not  have  to  take  the  inverse  function  as  we 
did  before  to  perform  the  regression  analysis.  Instead, 
we  take  the  square  difference  of  the  given  and  computed  t 
directly  from  this  expression.  Thus,  we  get 


'U 


.  ,  2  2 
a,+a.x, ,+a,x . . 

[Ty  -  exp {-10  1  2  *3  3  13  )]  , 


(29) 


for  i-th  absorption  band's  j-th  data  point,  where,  as 


before,  x^  is  given  by 


x  =  Cj  +  n  log(^J-)  +  m  logC^2-)  +  logU  .  (30) 

3  o  ij  3 

By  summing  this  individual  error  for  all  data  in  i-th  band, 
we  have  the  total  error  for  this  band  as 


E 


i 


E 

j  =  l 


(31) 


where  is  the  number  of  data  in  i-th  band. 

Again,  we  introduce  auxiliary  variables  u^,  i=l,2, 
...,NB  in  order  to  introduce  all  Cj,  i=l,2,...,NB  into  the 
Xjj  expression  Eq.(30).  By  this  we  get 


x 


ij 


NB  P,,  T 

l  .  uk,ijc’k  +  n  log(p  : )  +  m  + 

k= 1  ’  J  o  ij  J 


We  use  this  expression  for 


‘ij 


in  the  following  total  error 


(32) 

E 


E 


NB  J 
E  £  E 
i=l  j=l 


(33) 
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NB 

E 


i=  1 


Now,  we  are  ready  to  take  partial  derivatives  with 
respect  to  the  parameters  n,  m,  C j ,  C^B>  a^,  a2>  and 

a^  to  form  the  normal  equation  for  this  regression  problem. 
Theoretically  speaking,  we  can  evaluate  the  'best'  para¬ 
meter  values  by  solving  the  normal  equation.  But  obviously 
the  grand  total  Eq .  (33),  which  is  to  be  minimized,  is  not  a 
quadratic  function  of  the  unknown  parameters  and,  therefore, 
the  resulting  normal  equation  is  not  a  linear  function  of 
them.  Hence,  we  need  to  adopt  a  different  numerical 
method  for  the  evaluation  of  the  'optimal'  parameter  values. 

4,3  Nonlinear  Optimization  Method 

The  computational  technique  we  use  here  is  a  recursive 
technique  which  is  referred  to  as  the  conjugate  gradient 
method^  .  In  essence,  this  technique  improves  a  set  of 
guesses  of  the  parameter  values  recursively  by  locating  a 
new  set  of  guesses  which  yields  smaller  error.  For  a  given 
guess  (an,3n,  . . . ,yn)  of  the  minimizing  parameter  vector, 

at  which  the  error  is  minimized,  the  best  direction  of  the 
search  in  the  parameter  space  for  a  new  guess  is  first 
determined  using  up  to  second  order  derivatives  of  the 
error.  Then  the  one-dimensional  search  for  the  minimizing 
point  is  performed  along  this  direction  from  (otn,f3n,  ...» 
y  )  to  find  a  new  guess  (a  ,p  . ^  ) 
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which  yields  locally  the  smallest  error.  Now  this  procedure 

is  repeated  recursively  to  obtain  a  sequence  of  guesses 

until  the  gradients  become  less  than  a  small  positive 

number  which  is  chosen  a  priori. 

Actual  computation  was  done  by  utilizing  the  packaged 

1  A 

subroutine  FMCG  in  SSP  library  available  from  IBM  .  The 
necessary  gradients  are 

-2 ED  .  <5 f  .  , 
j  3 

-2 ED  6 f .  x  ,  , 

3  3  j 


3  J 
3  a . 


8  J 
3a , 


«  -2ED . 6f  x . 2 , 

3a3  3  j  3 

p 

=  -2EDj6fjCa2  +  2a3x.)log(?i) , 
I i  =  -2ZD.6fj(a2  +  2a3xj)log(^), 


(34) 


■  =  -2ZD.6f.(a_  +  2a,x  )u.  , 
9C ' .  332  3  j  i 

i 


NB  J 

where,  E  represents  E  E  and  D  and  6f  are  given  by 

i-1  j=l  J  J 


D 


j 


(35) 
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<$fj  =  (In  10)  10 


al  + 


l2xj  +  a3Xj 


fCx.), 


(36) 


and  f  (x)  is  given  by  Eq .  (28). 

Note  that  there  exists  a  linear  dependence  among  the 
gradients  which  is 


9  J 
l29a. 


+  2a 


9  J 
3  9a, 


NB 

E 

i=l 


(37) 


Therefore,  the  parameter  set  {n,  m,  Cj,  C^g,  a^,  a^, 

a^}  cannot  be  determined  uniquely.  As  it  was  explained  in 
the  previous  section,  this  is  due  to  the  arbitrariness  in 
the  positioning  of  the  standard  transmission  function. 
Hence,  the  spectral  parameter  Cj  is  again  set  to  be  zero, 
so  that  we  can  evaluate  unique  set  of  optimal  parameters. 
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V.  Comparison  of  the  Two  Methods 

5.1  Introduction 

Both  methods  can  evaluate  the  optimal  n,  m  and  C'(v) 
values  for  major  and  non-major  bands  and  also  a  standard 
transmission  function.  But  there  are  some  basic  differences 
which  are  discussed  in  the  sequel. 

5.2  Final  Products 

The  final  product  of  the  ADSET  code  is  a  piecewise 

analytical  standard  transmission  function  together  with 
the  band  model  parameters.  Each  analytical  piece  of  the 
standard  transmission  function  covers  only  one  of  the  pre¬ 
chosen  subintervals  of  (0,1)  transmittance  range.  On  the 
other  hand,  SIMMIN  produces  only  one  analytical  trans¬ 
mission  curve  for  the  entire  range.  Therefore,  ADSET  has 
more  flexibility  to  adjust  to  the  transmittance  curve 
variations.  This  feature  of  ADSET  can  be  very  valuable 
for  the  gases  with  non-standard  curves  of  growth. 

This  difference  is  amplified  when  the  number  of  the 
transmittance  sub-regions  used  in  ADSET  is  increased. 
However,  as  the  number  of  subregions  increase,  the  re¬ 
quirement  on  the  usable  data  becomes  severer  and  more 
spaces  i re  necessary  to  store  the  computed  results.  Hence, 
the  determination  of  the  number  of  subregions  should  be 
resorted  to  compromise. 
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5.3  Installation  of  the  Results  in  Lowtran 

The  final  products  of  two  codes  ADSET  and  SIMMIN  j 

were  installed  in  the  widely-used  Lowtran  code,  as  discussed 
in  Section  VII  of  this  report.  The  SIMMIN  results  require 
less  memory  space,  less  time  for  transmittance  computation 
and  simpler  coding  than  the  ADSET  result.  In  fact,  for  the 

SIMMIN  result,  all  that  have  to  be  stored  are  the  five 

band  model  parameters  n,  m,  a^,  a^,  and  a ^ ,  and  a  set  of 
spectral  parameters  C'(v^)  for  each  absorber.  Furthermore, 
the  computation  of  x  can  be  done  by  only  one  FORTRAN 
statement.  On  the  other  hand,  the  ADSET  result  requires 
the  storage  of  NCUT-1  of  a  ^  ,  a  ^  ,  and  a  values,  n  and  m 
and  a  set  of  C'(v)  values  for  each  absorber.  lhere  can  be 

a  large  difference  in  the  number  of  the  sets  of  a  ^  ,  a,, 

and  a.j  values  to  be  stored.  Moreover,  some  judging  state¬ 
ments  are  necessary  to  select  the  right  set  of  a^,  a,,,  and 
a  ^  [or  each  transmittance  computation. 

5.4  Data  Requirements 

The  ADSET  code  requires  the  cut  structured  data 
such  that  the  transmittance  of  each  data  ,ioint  mast  f  i[ 
in  one  of  prechosen  values.  But  the  SIMMIN  code  does  not 
impose  any  conditions  on  the  data  set. 

Some  considerations  on  the  requirement  of  equal 
transmittance  data  for  ADSET  are  due  here.  Even  if  the 
available  data  do  not  have  equal  transmittance  structure. 
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it  can  be  transformed  into  the  required  form  using  inter¬ 
polation/extrapolation.  This  constitutes  the  pre-processing 
of  the  raw  data.  Curves  of  growth  data  T  vs.  log  U  with 
T  values  not  necessarily  coinciding  with  the  prechosen 
values  can  be  locally  interpolated/extrapolated  using  an 
analytical  function.  This  procedure  is  indicated  in  Fig. 

3.  Again,  the  double  exponential  function  is  an  excellent 
choice  for  the  interpolation  function.  We  note  that  an 
almost  exact  technique  as  the  one  used  in  obtaining  a 
piecewise  analytical  transmission  function  can  be  used  for 
this  purpose.  In  fact,  only  a  minor  modification  of  the 
interpolation  subroutine  used  in  ADSET  can  accomplish  this 
task . 

5.3  Computation  Time 

The  numerical  methods  used  in  ADSET  and  SI MM IN  for 
solving  normal  equations  are  essentially  different.  The 
method  in  SIMMIN  is  a  recursive  algorithm  and  the  other  in 
ADSET  is  a  non-iterative  one.  Therefore,  the  computation 
time  for  ADSET  is  determined  by  the  size  of  the  data  set 
only,  whereas,  the  one  for  SIMMIN  depends  on  both  the 
actual  data  values  and  the  initial  guesses.  It  is  difficult 
to  estimate  the  computation  time  for  SIMMIN  due  to  this 
dependence.  One  way  of  controlling  the  time  is  to  limit 
the  number  of  iterations  performed.  This  feature  is  included 
in  the  packaged  subroutine  FMCG  which  is  used  for  actual 


26 


••PP-WWJI 


computations.  Actual  time  ranges  required  for  ADSET  and 
SIMMIN  computations  will  be  given  in  a  later  section. 


VI  . 


Lowtran  Capabilities  and  Functions 


6.1  Introduction 

The  Lowtran  code  consists  of  a  computer  model  for 
the  calculation  of  transmittance  through  atmospheres 
containing  absorbing  and  scattering  molecules  and  aerosols. 

The  models  used  in  the  code  were  for  the  most  part  developed 
in  1972^  but  later  editions  incorporated  computational 
changes  and  other  c ap abi 1 i t i e s ^  It  covers  the  spectral 

range  from  0.25  to  28.5  ym  at  intervals  of  5  cm  ^  with  a 
resolution  for  the  major  absorbers  of  20  cm  ^ .  The  trans¬ 
mittance  calculation  is  made  on  six  model  atmospheres  and 
two  haze  models  on  a  33-level  basis  for  altitude,  pressure, 
temperature  and  density  from  sea-level  to  100  km.  The  path 
of  the  transmission  is  considered  to  be  refracted  by  changes 
in  atmospheric  density,  a  fact  taken  into  account  in  an 
optional  subroutine.  In  its  present  form  the  Lowtran  code 
consists  of  a  single  main  program  that  inputs  the  path 
data  and  model  parameters,  computes  the  equivalent  absorber 
amount,  and  performs  the  transmittance  calculations.  The 
only  present  subroutines  are  associated  with  the  path,  and 
are  optional.  The  difficulties  of  understanding  and, 
especially,  updating  such  a  program  structure  are  considerable. 

The  principal  objective  of  this  effort  is  to  modify 
the  program  structure  of  the  Lowtran  code  in  accordance 
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with  the  following  criteria: 

1.  The  basic  functions,  calculations  and  print¬ 
outs  remain  nearly  identical  to  the  original. 

2.  The  basic  operations  involving  the  reading  of 
data,  the  calculation  of  the  equivalent  path 
and  the  transmittance  calculations  are  all 
separate,  independent  programs,  but  are  con¬ 
nected  as  subroutines  to  a  main  control  program. 

3.  The  structure  modification  is  performed  on  the 
latest  version  of  the  code,  i.e.,  Lowtran  4. 


As  an  exercise  in  the  use  of  the  modularized  version,  the 
present  authors  added  empirical  band  models  for  trans¬ 
mittance  through  the  trace  gases.  Also,  continuous 
functions  were  made  to  replace  their  transmission  tables 
for  the  principal  molecular  absorbing  species. 


r 
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6.2  Ceneral  Features 

In  this  section  an  effort  is  made  to  summarize  the 
basic  structure,  fundamental  calculations  and  models  used 
in  the  Lowtran  code  for  estimating  atmospheric  attenuation 
by  gases  and  aerosols.  Reference  is  specifically  made  to 
the  latest  fourth  version,  although  at  the  present  time 
the  authors  are  aware  of  a  recent  effort  by  AFGL  on  a 
fifth  version.  From  the  authors'  evaluation  of  their 
recent  efforts,  it  appears  that  the  modularization  presented 
here  may  be  incorporated  in  their  latest  version.  For 
instance,  the  latter  is  known  to  have  a  single  separate 
subroutine  for  the  emission  and  transmission  calculations. 
The  major  contribution  of  the  work  presented  here  lies 
in  the  separation  of  that  emission  and  transmission  loop 
into  a  subroutine  for  model  selection,  a  subroutine  for 
the  equivalent  path  and  individual  subroutines  for  all  of 
the  attenuation  models  in  the  code. 

The  Lowtran  code  is  designed  for  the  specific 
purpose  of  calculating  at  low  resolutions  either  atmospheric 
radiance  or  transmittance  between  any  two  locations  in 
the  Earth's  atmosphere  at  frequencies  ranging  from  the 
ultraviolet  (UV)  to  the  infrared  (IR).  This  is  accom¬ 
plished  through  the  use  of  band  models  accounting  for 
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resonant  gaseous  absorption  (e . g .  H^O  vapor,  0^  ,  HNO^ 
vapor  and  the  uniformly-mixed  gases),  resonant  aerosol 
absorption,  non-resonant  gaseous  absorption  (e . g .  ^  and 

H^O  vapor  continua)  and  scattering  by  molecules  and  aerosols. 
The  spectral  intervals  over  which  the  band  models  are 
provided  vary  from  5  cm  ^  to  500  cm  ^ ,  as  shown  in  Table  1. 

It  should  be  pointed  out  that  the  spectral  resolution  is 
generally  much  lower  than  the  interval  over  which  they 
are  defined.  For  instance,  the  models  for  the  principal 
absorbers  are  given  at  5  cm  ^  intervals,  while  their 
spectral  resolution  is  20  cm  ^ .  The  spectral  resolutions 
for  the  remaining  models  is  not  specified  anywhere  in  the 
available  literature  on  the  code.  In  this  table  the 
spectral  definition  of  the  models  for  aerosol  absorption, 
and  for  aerosol  and  molecular  scattering  are  not  shown 
because  they  are  spectrally  continuous. 

The  spectral  regions  over  which  the  attenuation 
models  are  effective  are  summarized  in  Table  2.  It  may 
be  seen  in  this  table  that  over  some  regions  only  a  few 
species  attenuate  and,  therefore,  a  transmittance  of  unity 
may  be  specified  in  the  calculation  of  the  total  trans¬ 
mittance.  This  table  forms  the  basis  for  the  model 
selection  subroutine  introduced  in  the  modularized  version 
for  the  purpose  of  simplifying  the  code  structure. 

In  the  discussion  that  follows,  the  individual 


Frequency  interval  of  the  attenuation  band  models 
in  the  Lowtran  code.  The  models  for  aerosol 
absorption  and  aerosol  and  molecular  scattering 
are  spectrally  continuous  and,  therefore,  not 
shown . 


ATTENUATING  WAVENUMBER  SCALE  (cm1) 

SPECIE  1000  2000  3000  4000  5C 


LOWTRAN  are  effective 
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attenuation  models  are  grouped  together  in  certain  classes 
and  are  briefly  discussed.  Generally  speaking,  the  dis¬ 
cussion  is  restricted  to  the  extent  of  illustrating  the 
function  and  parameters  which  had  to  be  identified  in 
Lowtran  for  the  modularization  purpose  that  followed.  An 
exception  is  made  in  the  case  of  the  major  molecular 
absorption  models  (i.e.  1^0  vapor,  infrared  0^  and  the 

uniformly-mixed  gases)  because  they  are  replaced  with 
continuous  functions  in  the  modularized  version.  For  a 
comprehensive  discussion  on  the  theory  of  all  of  the 
original  models  the  reader  is  encouraged  to  study  the 
series  of  AFGL  reports^  ^  on  the  code,  as  well  as  the 
references  therein. 

6.3  Resonant  Molecular  Absorption  Models 

Molecular  resonant  absorption  is  modeled  in  the 
code  for  vapor,  infrared  0^  ,  the  uniformly-mixed  gases, 

and  HNO^  vapor.  Different  approaches  are  used  for  the 
first  three  listed  as  compared  with  the  approaches  used 
in  connection  with  0^  in  the  visible  and  ultraviolet 
regions  and  with  HNO^  vapor  in  the  infrared. 

The  models  used  to  account  for  gaseous  absorption 
by  the  molecules  of  1^0  vapor,  infrared  0^ ,  and  the 
uniformly-mixed  gases  are  based  on  Eq.(8),  namely 


T  =  f  {  X  }  . 


(8) 
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The  developers  of  Lowtran  obtained  the  parameters  n,  m, 
the  function  f  and  the  spectral  constant  C'  at  5  cm  1 
intervals  using  experimental  and  calculated  transmittance 
data  of  20  cm  ^  resolution.  Table  3  shows  the  values  of 
the  parameters,  as  well  as,  the  equations  for  the  calcu¬ 
lation  of  th  ’  absorber  amount.  The  spectral  constant 
C*  over  the  entire  spectrum  of  definition  may  be  found  as 
part  of  the  data  input  presented  in  the  Appendix.  The 
transmission  model  for  the  uniformly-mixed  gases  was 
obtained  by  combining  the  data  for  all  of  these  gases  in 
the  proportions  listed  in  Table  4.  It  should  be  pointed 
out  that  the  temperature  and  pressure  exponents  used  in 
Lowtran  for  the  major  absorbers  and  listed  in  Table  3  are 
not  the  same  as  the  ones  developed  from  the  original 
transmission  data.  This  inconsistency  was  introduced 
during  the  digitizing  of  the  curves  for  inclusion  in  the 

computer  code,  in  order  to  account  more  accurately  for 

2 

the  temperature  dependence  . 

The  method  used  for  modeling  HNO^  vapor  and  the 
visible  and  ultraviolet  0^  is  similar  to  the  one  described 
above  for  the  major  absorbers,  except  that  the  function 
f  was  specified  a  priori  to  be  an  exponential.  Thus, 


where  for  HNO^ 


x  =  exp (-CW) , 


(38) 
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o 


U  =  M  Z  x  105, 


and  for  0^ 


W  =  U  =  46.667  pz. 


(39) 

(40) 


(41) 


In  Eq.  (40)  M  is  the  mixing  ratio  profile  as  tabulated 
in  the  Appendix  together  with  the  C '  s,  and  p  is  the  absorber 
density . 

The  last  of  the  molecular  absorption  models  is  the 
one  for  the  resonant  absorption  by  atmospheric  aerosols. 

The  exponential  function  i ..  Eq  .  (38)  is  assumed 

T  =  exp (-CW) , 


where 


W  =  U  =  3.5336 


x  10_6NZ, 


(42) 


and  N  is  the  vertical  distribution  of  the  number  of  haze 
particles.  Tabulations  are  provided  of  distributions  for 
5  Km  and  23  km  visibility,  as  listed  in  the  Appendix.  Other 
visibilities  are  treated  in  the  code  itself  through  linear 
interpolat ion . 


6.4  Non-Resonant  Molecular  Absorption  Models 

Non-resonant  gaseous  molecular  absorption  is  re¬ 
presented  by  the  N£  and  H^O  vapor  continuums.  The  same 
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modeling  approach  is  used  for  N.;  as  for  resonant  molecular 
absorption,  that  is 

T  =  exp (-CW) , 

where 

p  2  T  1.5 

W  =  qb)  (^)  U  (43) 

o 

U  =  0.8  Z.  (44) 

For  the  H^O  vapor  continuum  an  exponential  function  is 
also  used,  but  with  a  more  elaborate  exponent.  Thus, 


T  =  exp ( -y ) , 


(45) 


where 


C 

y  =  C  [P  +  ~  (P  “  P  ) ]  U.  (46) 

s  w  C  w 

s 

Here,  P^  is  the  partial  pressure  of  water  and  Cg  and 
are  the  self-broadening  and  nitrogen-broadening  spectral 
constants.  The  values  of  these  spectral  constants  depend 
on  the  spectral  region  where  the  continuum  is  effective. 
In  the  8  to  14  pm  region 

C  =  C  exp  [  6 . 08  (•—  -  1)]  ,  (47) 

so  1 


and 

C 

-jr-  =  0.002,  (48) 

s 

while  in  the  3.5  to  4.2  ym  region 


r 


39 


Cs  =  Co  exp  [  4 . 56  ^  -  1)]  , 


and 


C 

=  0.120. 

S 


In  these  equations  the  value  of  is  given  by 
CQ  =  4.18  +  5578  exp(-7.87  x  10_3v). 


'■  i 


I 


(49) 


(50) 


(51) 
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ATTENUATING 

SPECTRAL 

PRESSURE 

TEMPERATURE 

ABSORBER 

SPECIE 

REGION 

EXPONENT 

EXPONENT 

AMOUNT 

(cm-1) 

n 

m 

U 

Ho0  Vapor 

350-  9,195 

i 

L. 

9,875-12,795 

0.90 

0.45 

0.1  pZ 

13,400-14,520 

Uni f  o  rmly- 

i 

Mixed  Gases 

500-  8,070 

1.75 

1.375 

Z 

12,950-13,245 

Infrared 

°3 

575-  3,270 

0 .40 

0.20 

46.667  pZ 

N, 

Cont Inuum 

2,080-  2,740 

2 .00 

1 . 50 

0.8  Z 

Aero  sol 

Absorpt ion 

333-50,000 

0.00 

0.00 

3.5336  x  10~6  NZ 

Aerosol 

Scattering 

333-50,000 

0 . 00 

0.00 

3.5336  x  10~A  NZ 

Molecular 

-20 

9.87  x  10  U  Z 

Scattering 

3,000-50,000 

1 . 00 

1 . 00 

hno3 

850-  920 

Vapor 

1,275-  1,350 

1 . 00 

1 . 00 

1  x  105  MZ 

1,675-  1,735 

Visible  and 

13,000-24,000 

Ultraviolet  0^ 

2  7 ,500-50,000 

0.00 

0.00 

46.667  pZ 

Table  3.  Absorber  parameters  in  I.owtran  for  the  attenuation  models, 

where  p  is  the  density,  Z  the  range  and  M  the  mixing  ratio 
The  H_0  continuum  model  is  excluded  because  of  its  dif¬ 
ferent  functional  form. 
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GAS 

MOLECULAR 

WEIGHT 

PARTS  PER  MILLION  : 
VOLUME  (ppm) 

co2 

44 

330.0 

n2o 

44 

0.28 

CO 

28 

0.075 

ch4 

16 

1 .60 

°2 

32 

2.095  x  105 

Table  4.  Concentrations  of  the  uniformly-mixed  gases 

used  in  the  combined  model. 


42 


6.5  Scattering  Models 


In  order  to  account  for  atmospheric  scattering 
exponential  functions  were  used  again.  For  scattering 
by  molecules  the  model  is  defined  as  in  Eq .  (38) 

X  =  exp(-CW) 


where 

T 

W  =  (f-)(^)  U  (52) 

o 

U  =  9.87  x  lu-20  Z  (53) 

C  =  V4  (54) 


For  aerosol  scattering  the  argument  of  the  assumed  ex¬ 
ponential  function  is 


W 


U 


3.5336  x  10 


NZ 


(55) 


VII  . 
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Modularization  of  Lowtran  Including  the  Trace  Gases 

7.1  Introduction 

Considering  the  generality  and  broadness  in  scope 
of  this  code  it  is  not  surprising  that  the  program  structure 
shows  in  its  present  form  great  complexity.  Although  the 
program  user  is  not  normally  interested  in  aspects  of  the 
code  other  than  the  input  and  output,  there  are  many  cases 
where  a  basic  understanding  helps  in  specific  applications. 
Situations  are  likely  to  occur,  for  instance,  where  a 
replacement  of  one  of  the  several  attenuation  models  is 
highly  desirable.  To  assist  in  the  implementation  of 
model  additions  or  changes  as  well  as  in  the  extension  to 
other  spectral  regions  and  media,  the  concept  of  the 
modularized  version  was  conceived.  This  version'*'"’  was 
designed  to  represent  exactly  the  same  calculations  as  the 
original,  except  for  the  simplification  of  the  program 
structure  into  modules  or  subroutines.  However,  upon  the 
termination  of  that  task  the  authors  proceeded  to  add 
models  for  the  trace  gases,  as  developed  during  the  pre¬ 
sent  scientific  effort. 

7.2  Structure  of  Modularized  Version 

The  basic  design  used  was  that  of  a  main  program 


which  reads  input  data,  computes  total  transmittance  and 
radiance  and  generates  outputs,  and  a  series  of  subroutines 


which  select  individual  models  and  compute  individual  trans 
mittances  and  absorber  amounts.  This  is  shown  in  Fig.  A. 
The  main  operational  flow  chart  follows  in  Fig.  5.  Ex¬ 
cluding  the  four  subroutines  for  the  trace  gases,  the 
modularized  version  breaks  down  the  original  into  one 
program  with  11  subroutines.  The  flow  chart  for  subroutine 
ABSORB  is  shown  in  Fig.  6.  This  subroutine  computes  the 
equivalent  absorber  amount  for  all  of  the  attenuation 
models  according  to  Eq.(^),  which  in  terms  of  the  meteo¬ 
rological  variables  becomes 

p  r  y  v  o  T  m 

w  (uzj)  dU  (56) 

Figure  7  gives  details  of  the  Transmittance/Radiance  Loop 
of  program  Main.  It  is  worth  noting  that  the  modularized 
version  of  Lowtran  being  done  by  AFGL  separates  this  loop 
into  a  subprogram.  The  modularization  discussed  in  this 
text  leaves  the  loop  as  part  of  the  main  program,  but 
extracts  individual  subroutines  for  the  calculation  of  the 
equivalent  absorber  amount,  the  frequency  selection,  and 
the  attenuation  models. 

The  flow  chart  for  FREOSL  subroutine  is  shown  in 
Fig.  8.  This  subroutine  is  designed  to  simplify  the  pro¬ 
cess  of  arriving  at  the  individual  models  effective  at 
the  frequency  of  interest.  It  should  also  assist  the 


ABSORB : 

Computes  absorber 
amounts 


MAIN  : 

Reads  basic 
model  data , 

v, 

FREQSL : 
Selects 

models 

1 

c  ompu  tes  total 

r 

effective  at  given 

transmittance  & 

frequency 

J 

radiance  ,  prints 
outputs 

> 

1 

Thirteen  model 
subroutines 

ig.  4 


Conceptual  flow  chart  of  modularized  Lowtran. 
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Program  Main 

Read  in  basic  data  model 
atmospheres  transmission 
functions  spectral  data 


Entry  Point  2 


Read  in  and  print  problem 
parameters  MODEL,  IHAZE , 

I TYPE ,  etc. 

.  ~T . . . . 

If  IEMISS  =  1  print  108 
If  IEMISS  =  0  print  109 


Entry  Point  3 


If  IHAZE  =  7 


Read  in  and  print 
aerosol  model 


Entry  Point  5 


r 

_ 

ABSORB 

bllbf OUClllti  AbbURo 
Compute  equivalent 
absorber  amount 

Entry  Point  6 


Print  equivalent  path  quan¬ 
tities  for  each  absorber 
W(K) ,  K  =  1 ,  KMAX 


X 


Transmittance/Radiance 

LOOP 

- - T - 


Subroutine  FREQSL 
selects  13  models 


Read  IXY 


If  IXY  =  0 


> 


END 


Read  VI,  V2 ,  DV 
for  same  path 


If  IXY  =  1 


If  IXY  =  2  ^>- 

X 


li¬ 


st  art  new  problem 


If  IXY  =  3  ^>- 

I 

I> 


If  MODEL 


> 


If  IXY 


Read  MODEL,  IHAZE, 
etc.  for  same  path 


Fig.  5. 


General  flow  chart  for  Modularized  Lowtran  4 


Subroutine  ABSORB 


If  IXY 


If  IXY 


Read  in  and  print  Hi,  H2, 
ANGLE  RANGE  BETA,  VIS,  VI, 
V2,  DV 


If  MODEL 


If  MODEL 


If  BETA  =  0  CALL 
ANGLE  Subroutine 


Read  in  and 
Print  met.  data 


Read  radiosonde 
data  ML  levels 


Subroutine 
ANGIE  Calculate 
ANGLE,  given 
HI,  H2,  BETA 


ANGLE 


Set  up  "Horizontal  Profiles"  EH(K,I) 
i.e.  Equivalent  Path  Quantities  per 
km  path  length  at  each  level(I)  and 
for  each  absorber  (K) 


Determine  EH  Quan¬ 
tities  at  Altitudes 
HI  and  H2 


HI.  H2 


EH  at  X 


Set  Total  Slant  Path  Absorber 
Amounts  for  Each  Absorber  to 
Zero,  i.e.  W(K)=0  for  K=l,  KMAX 


Subroutine 


POINT 


Fig.  6,  Flow  chart  for  subroutine 

ABSORB,  computing  equivalent 
path  quantities 


©  © 


If  ITYPE-1 


If  ANGLE>90 


Horizontal  path 
HI  and  RANGE 
given 


W«RANGE*EH 


Determine  HMIN 
3g  minimum  height 

- -  of  path  also 

refractive  index 
above  HMIN 


Slant  path  to 
Space:  Hi  and 
ANGLE  given 

27 

Calculate  Slant  Path 
Quantities 

and  0  for  each  layer 
path 


L 


32  Load  WLAY  Matrix 


58  I 

Return 


Slant  Path  from 
HI  to  H2  at 
ANGLE 


Calculate  Slant 
5  path  quantities 
and  0 

for  each  layer 


H2=HMIN 


-HI 

w=  J 

EH 

H2 

50  Load  WLAY  Matrix 
- 1 - - 


If  path  does  not  pass 
through  HMIN  or  if 
H2  -  HMIN 


- , - 

H2 

W=W+2  /  EH 

HMIN 

H2 

W=2  /  EH 
HMIN 


Reset  ANGLE  to  be  (180-ANGLE  + 
refraction  contribution)  and 
proceed  as  for  an  upward 
looking  case 


Fig.  6. 


(cont'd) 
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Input  equivalent  absorber  W(K) 
for  required  atmospheric  path 
Print  W(K){  K-l,  KMAX 


If  Vl  and/or  V2  are  set  outside 
of  the  range  of  the  prograij  that 
are  reset,  e.g.  Vl=350  cm  freq¬ 
uency  V  initially  set:  V=V1-DV 


Subroutine  PATH  compute 
matrix  WPATH 


if  IEMISS=U 
1KMAX=15 


Do  17  IK=1 ,  IKMAX 


CALL 

FREQSL 


SUBROUTINE  FREQSL 
determine  effective 
transmittance  models 
for  a  given  frequency  V 


Thirteen  Transmittance 
model  subroutines 
compute  TX(l)->TX(15) 
for  a  given  frequency  V 


Tx(9) =SUM4  +  SUM5  +  SUM6  +  SUM7 
Tx(9)**EXP{-Tx(9) } 

Tx(9)»Tx(l) *Tx(2)*Tx(3)*Tx(9>* 
Tx  ( 12 )  *Tx  ( 13)  *Tx  ( 14)  *Tx  ( 15 ) 
Tx(9)  is  set  equal  to  total 
transmittance 


Tx(10) =1-Tx(10) 
SUHA-ETx(9)*DV 


Output  total  radiance 


Print  V.A.Tx(l)-*Tx(7) tTx(10) 


If  V.GE.V 


V-V+DV  I - 1  Print  V, A ,Tx(ll)~Tx(15) ,Tx(9) 


Read  IXY 


See  Fig.  6  . 


Fig.  7.  Flow  chart  for  transmittance/radiance  loop. 
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If  (I.GE.01  and 


I.  I.E.  13, 


If  (I.GE.19  and  I.LE.30) 


If  (I.GE.3J  and  I. I.E. 45), 


If  (I. OF..  56  and  I.LE.5A) 


If  (I.CE.55  an  d  I.  I.E.  60) 


T'VhOr 


If  (I.G”.61  and  I. IE. 64) 


Fig.  8  .  flow  chari 


Subroutine  I.IJAP 
P=C+LofW(l) 

Tx( 1) =EXP [-10**(-1, ] 4619  + 
0.55013*p) ] 


Go  to  40 


CALL 

LUAP 

EKUL 

SUSEJ 


FUSEJ 


Subroutine  DTVAD 
Uniformly  mixed  gases 
P-C2+LogW(2> 

Txf2)=EXP[-10**(-l. 14619  + 
0.5501  3*p)  ] 


Go  to  40 


CAM. 

LUAP 

DTVAI) 

EVEN’S 


G 

o  to  41 

CALL 

Subroutine  L'VETS 

LEAP 

Or  one 

!)' VA0 

?=-C3+LogT(3) 

LVFTS 

Tx(3)«l/[l+F.XF(-3. 08019  + 

i'Kl'l. 

2. ll]27*p) 1 

Go  to  -tO 


.r.vp 

n>  cad 

i.Vl.TS 

..Kl'L 


.'.Kl' 

IMJ 


Subroutine  KPAM 
15,  continuum 
Tx(4)«C4*V(4) 
S'.2%«T*(4) 
T:,(4)«KVr(-Tx(4)  } 


, >  to  40 


subroutine  'TJiQSL  and  thirteen  transmittance 
nod’)  subroutines. 


5 


Go  to  40 


ig.  8< 


Continued 


■ 
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Fig.  8. 


Conti  nufed 


If 

(T.GE.278 

and 

I. LE, 282)  \ 

or 

(I.GE, 326 

and 

I.LE.346)  \ 

or 

(I.GE. 480 

and 

I.LF.498)  / 

or 

(I.GE. 511 

and 

I.LE.531)  / 

If  (I.GE. 283  and  I.LF.325) 


CALI. 

LUAP 

DI.VAD 

EVETS 

NHOJ 

EKUL 


CALL 

LUAP 

DIVAD 

EVETS 

NHOJ 

EKUL 

ARZE  i 


Subroutine  Sl’SEJ 
SO, 

PaC1 2+LogW(12) 
Tx(12)-EXP[-10**(51  + 
52*p) ] 


Go  to  40 


Subroutine  ARZE 
NO 

P=C1 3+Lo  gW( 1 3) 

Tx(  13)  =EXP  f- 1.0**(  FN1  + 
FN2*p) ] 


Go  to  40 


If  (I.GE.347  and  I.LE.420) 
or  (I.GE.439  and  I.LE.479) 


CALI: 

LUAP 

DIVAD 

EVF.TS 

KRAM 

NHOJ 

EKUL 


Subroutine.  St’TIT 

nh3 

P®C14+LoojW  (14) 
Tx(14)=SXPi-10**CFHl  + 
FH2*p) ] 


Go  t o  40 


If  (I.GE.421  and  I.I.E.438) 


CALL 

LUAP 

DIVAD 

EVETS 

KRAM 

NHOJ 

EKin 

SUSET 


Subroutine  ROJ 
N°2 

P=C15+LogW(15) 

TxO  5)  “EXP  [-10**(0l  + 

02*p)] 


user  who  desires  to  replace  or  add  models  to  the  program 


and  it  should  reduce  the  overall  computational  time.  This 
subroutine  is  based  on  Tables  2  and  8, 

7.3  Models  for  H?0  Vapor,  Infrared  0„  and  the  Uniformly- 
Mixed  Gases 

As  indicated  abovd,  all  the  attenuation  models  were 

extracted  from  the  main  program  and  placed  into  subroutines. 

The  models  were  left  basically  in  the  same  structural  form 

except  for  the  models  for  HNO^  and  vapor,  infrared  0^ 

and  the  un i f o rml y -mixed  gases.  The  change  in  the  first 

was  to  arrange  it  along  the  same  form  as  in  the  other 

models  originally  available  in  Lowtran.  That  is,  the 

spectral  parameters  were  extracted  from  the  subroutine 

and  read  at  the  beginning  of  program  MAIN.  The  changes 

in  the  latter  three  gases  (i.e.  H^O  vapor,  0^  in  the 

infrared  and  the  uniformly-mixed  gases)  were  based  on  a 

2 

previous  work  by  Pierluissi  et  al.  on  the  representation 
of  the  tabulated  transmission  functions  by  analytical 
functions.  The  other  principal  change  consisted  of  adding 
models  for  the  trace  gases  S  0  ^  ,  NO,  NO  ^  and  N  H  ^  • 

To  arrive  at  the  analytical  function  for  modeling 
1^0  vapor  and  the  uniformly-mixed  gases  the  double  expo¬ 
nential  expression 

T  =  exp (-10ao+alX)  (57) 

where  x  is  as  in  Eq.  7  and  a  and  a.  are  absorber 

^  o  1 


i 


constants,  was  curve-fitted  to  the  134  values  of  x  and  x 

tabulated  in  Lowtran.  The  values  found  are  a  =  -1.14619 

o 

and  a^  =  0.55013,  and  it  reproduced  the  tabulated  trans¬ 
mittance  with  a  standard  deviation  of  0.005.  For  0^  the 
function  adopted  is  given  by 


T 


_1 _ 

a  +  a  ,  x 
o  1 


(58) 


where  a^  =  -3.08019,  a^  =  2.11127,  and  the  tabulated  data 
is  reproduced  with  a  standard  deviation  of  0.007.  Note 
in  each  one  of  these  functions  that  the  134  tabulated 
values  are  replaced  with  two  and,  hence,  their  adoption 
reduces  the  computer  storage  requirements.  Also,  they 
inherently  offer  exponential  interpolation  while  with 
the  present  tabulation  linear  interpolation  is  being 
used.  Finally,  tii  ore  is  no  need  for  the  small  optical 
thickness  (i.e.  0.999  ■'  T  <  1)  correction  inserted  in 

Lowtran  4,  as  required  by  its  radiance  calculational 
s  cheme . 


7.4  Models  for  Trace  Cases  SO^,  NO,  NO^,  and  N  H  ^ 

Absorption  by  the  trace  gases  was  incorporated  in 
Lowtran  using  a  somewhat  similar  procedure.  Empirical 
transmission  functions  were  first  obtained  from  a  com¬ 
puterized  procedure  which  replaced  the  classical  manual 
graphical  techniques.  The  procedure  is  explained  in 

Chapter  III  of  this  report  and  has  been  proposed  to  the 

1  1 


scientific  community 
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Instead  of  either  representing  the  transmission  function 
by  a  table  or  by  a  single  function,  it  was  divided  into 
nine  segments  for  each  absorber.  The  individual  curve 
segments  are  summarized  in  Table  5,  each  one  being  repre¬ 
sented  by  the  function 

a  +  ax 

T  =  exp (-10  u  1  )  (59) 

For  each  absorber  x  is  computed  with  Eqs.  (8)  through  (11) 
and  the  relation 

U  =  0.772  x  10  4  ppm  p  Z  (60) 

a 

where  ppm  is  the  parts  per  million  by  volume,  p  is  the 

d. 

3 

air  density  in  gm/m  and  Z  is  the  range  in  kilometers. 

Table  6  lists  the  ppm  and  temperature  and  pressure  exponents 
used  in  the  modularized  code  for  the  individual  trace 
gases.  The  ppm  values  are  read  as  input  through  a  separate 
card  which  may  be  easily  changed  according  to  the  needs  of 
the  user.  The  constants  C'  are  tabulated  in  Table  7.  The 
spectral  coverage  for  each  gas  is  depicted  in  Table  8. 

The  models  are  for  a  resolution  of  20  cm  ^  and  are  defined 
at  5  cm  ^  through  their  spectral  regions  of  effectiveness. 
Their  mean  standard  deviation  in  fitting  the  original  line- 
by-line  data  is  about  0.008.  Figure  9  depicts  the  trans¬ 
mission  functions  for  the  four  trace  gases  considered. 


WAVENUMBER 

( cm--*- ) 

C 

I 


WAVENUMBER 

(cm-l ) 

C 

-2.987 

-2.330 


-1.370 
-1.041 
-0. 795 
-0.613 
-0.469 
-0.346 
-0. 233 
-0.126 
-0.037 
0.0 
-0 . 008 
-0.052 
-0 . 102 
-0.102 
-0.044 
0. 013 
0.039 
0.014 
-0 .056 
-0. 14  L 


1070 
1075 
1080 
1085 
1090 
1095 
1100 
1105 
1110 
1115 
1120 
.25 
1130 
1135 
1140 
1145 
1150 
55 
60 
65 
1170 


-1.080 
-0.926 
-0.787 
-0 . 661 
-0. 544 
-0.434 
-0.329 
-0 .230 
-0.139 
-0.073 
-0.047 
-0.057 
-0.083 
-0 . 098 
-0.071 
-0.020 
0.014 
0.011 
-0.040 
-0.123 


WAVENUMBER 
(cm-1 ) 


1320 
.325 
1330 
1335 
1340 
.345 
1350 
355 
.360 
365 
1370 
137  5 
1380 
1385 
1390 
1395 
1400 
14  0  5 
1410 
2450 
2455 
2460 
2465 


1 


-0.494 
0 .139 
0.613 
0.899 
1 . 043 
1 . 090 
1.097 
1.104 
1 . 093 
1.118 
1.088 
0.926 
0.534 
-0.067 
-0.804 
-0 .768 
-1 .687 


2.855 


■2  . 

•1 .528 
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WAVENUMBER 

(cm"-*-) 

m 

WAVENUMBER 

(cm-1) 

n 

WAVENUMBER 
(cm- 1 ) 

H 

555 

-0 .221 

1185 

-0.213 

2470 

-1.076 

560 

-0.294 

1190 

-0 . 301 

2475 

-0.805 

565 

-0.366 

1195 

-0. 388 

2480 

-0.647 

5/0 

-0.442 

1200 

-0.481 

2485 

-0.571 

575 

-0.529 

1205 

-0.586 

2490 

-0.549 

580 

-0.635 

1210 

-0.707 

2495 

-0.539 

585 

-0.766 

1215 

-0.843 

2500 

-0.536 

590 

-0.934 

1220 

-0. 996 

2505 

-0.517 

595 

-1.157 

1225 

-1.165 

2510 

-0. 528 

600 

-1.457 

1230 

-1.351 

2515 

-0.691 

605 

-1 .862 

1235 

-1 .554 

2520 

-1.073 

610 

-2.420 

1240 

-1.777 

2525 

-1.673 

615 

-3 .094 

1245 

-2.033 

2530 

-2.414 

1055 

-2 .604 

1250 

-2 .369 

2535 

-2.207 

1060 

-2.156 

1310 

-3 .010 

1065 

-1.884 

1315 

-2.080 

Table  7a 


(Cont inued ) 


820  -0.519  I  1895  -0.151  1970  -2.806 


c ' 

WAVENUMBER 

-0.844 

800 

-0.760 

805 

-0.676 

810 

-0.60b 

815 

-0.543 

820 

-0.496 

825 

-0. 450 

830 

-0.414 

835 

-0 . 383 

840 

-0.326 

845 

-0.289 

850 

-0.217 

855 

-0  .  140 

860 

-0 . 097 

865 

-0.034 

870 

-0. 031 

875 

-0 . 082 

880 

-0.139 

1,540 

-0.21b 

1  ,545 

-0.249 

1,550 

-0.207 

1,555 

-0.117 

1 ,560 

-0.255 

-0.286 

-0.315 

-0.334 

-0.352 

-0.36b 

-0.396 

-0.423 

-0.459 

-0.498 

-0.541 

-0.586 

-0.630 

-0.676 

-0.720 

-0.766 

-0.809 

-2.428 

-1.494 

-0.647 

0.122 

0.756 


1 ,600 
1,605 
1,610 
1,615 
1,620 
1,625 
1 ,630 
1,635 
1 ,640 
1,645 
1  ,650 
1 ,655 
1 , 660 
1,665 
1,670 
2,840 
2,845 
2,850 
2,855 
2,860 
2 ,865 
2,870 


2.576 
2  .  350 


-1 

220 

-0 

644 

-0 

253 

0 

052 

0 

326 

0 

574 

0 

792 

fficient  C'(v)  for  NO2  ■ 


6 


C  ' 

WAVENUMBER 

| 

WAVENUMBER 

C  ' 

765 

-0.047 

1,565 

1.230 

2,875 

0.978 

770 

0 . 000 

1,570 

1.568 

2,880 

1.122 

775 

0.009 

1 .855 

2,885 

1  .  216 

780 

-0.046 

■SB 

2 .104 

2,890 

1.252 

785 

-0.100 

■sb 

2 .310 

2,895 

1 . 249 

790 

-0.148 

■SB 

2.469 

795 

-0.214 

■ 

2.573 

Table  7c. 
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WAVENUMBER 

C  ' 

WAVENUMBER 

' 

wavenumber 

C  ' 

690 

875 

-1 

.124 

1,060 

-0 

.  589 

695 

-2 

880 

-1 

.155 

1,065 

-0 

.565 

700 

-2 

.290 

885 

-1 

.161 

1,070 

-0 

.537 

705 

-2 

.128 

890 

-1 

.143 

1,075 

-0 

.  510 

710 

-1 

.980 

895 

-1 

.139 

1 ,080 

-0 

.  512 

715 

_  2 

.225 

900 

-1 

.117 

1 , 085 

-0 

.  528 

720 

-1 

.823 

905 

-1 

.107 

1,090 

-0 

.575 

725 

-1 

.744 

910 

-0 

.844 

1 , 095 

-0 

.625 

730 

-1 

.674 

915 

-0, 

.  558 

1 ,100 

-0. 

.  668 

735 

-1 

.577 

920 

-0. 

.  238 

1,105 

-0  , 

.  694 

740 

-1  . 

.481 

925 

-0, 

.042 

1,110 

-0. 

.717 

745 

-1  . 

.  372 

930 

-0. 

.  002 

1,115 

-0  , 

.740 

750 

-1  . 

.284 

935 

-0  . 

.  157 

1,120 

-0. 

.774 

755 

-1  . 

.207 

940 

-0. 

,436 

1 , 125 

-0. 

,834 

760 

-1  . 

,128 

945 

-0. 

,610 

1,130 

-0. 

,905 

765 

-1  . 

,061 

950 

-0  . 

548 

1,135 

-0. 

977 

770 

-1  . 

004 

955 

-0  . 

352 

1,140 

-1 . 

042 

775 

-0. 

947 

960 

-0  . 

139 

1,145 

-1 . 

133 

780 

-0  . 

886 

965 

-0  . 

09  5 

1,150 

-1 . 

219 

785 

-0  . 

876 

970 

-0. 

365 

1,155 

-1 . 

301 

790 

-0. 

872 

975 

-0. 

729 

1,160 

-1 . 

383 

795 

-0  . 

869 

980 

-1 . 

048 

1,165 

-1 . 

488 

800 

-0. 

872 

985  j 

-1 . 

275 

1,17  0 

-1 . 

594 

Table  7  d . 

The 

spectral  coefficient  C'(v)  for  NH^. 

WAVENUMBER 

C  ' 

WAVENUMBER 

C  ' 

WAVENUMBER 

mm 

805 

-0. 848 

990 

BB 

1,175 

810 

-0.811 

995 

Hi 

1,180 

-1.796 

815 

-0.772 

1,000 

-1 . 053 

1,185 

-1.873 

820 

-0.773 

1,005 

-0. 963 

1,190 

-1 .936 

825 

-0.793 

1,010 

-0. 920 

1,195 

-1 . 991 

830 

-0.825 

1,015 

-0.944 

1 , 200 

-2.080 

835 

-0.869 

1 ,020 

-0.889 

1,205 

-2  .  183 

840 

-0.894 

1,025 

-0.829 

1,210 

-2.292 

845 

-0.890 

1,030 

-0.736 

1 ,215 

-2.404 

850 

-0.873 

1 , 035 

-0.644 

1,220 

-2.529 

855 

-0.868 

1,040 

-0.596 

1,225 

-2.639 

860 

-0.907 

1,045 

-0.569 

1  ,230 

-2.732 

865 

-0.965 

1,050 

-0.572 

870 

-1.045 

1,055 

-0.590 

Table  7  d  . 
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VIII.  Calculations  and  Results 
8.1  Introduction 

The  procedure  for  the  use  of  the  Modularized  Lowtran 
in  calculations  is  identical  to  that  of  the  original  and, 
hence,  deserves  no  further  explanation.  There  are  some 
input  and  output  alterations  that  deserve  some  explanatory 
remarks.  Changes  in  the  input  format  include: 

1.  Reading  of  the  spectral  constants  for  all  band 
models  at  the  beginning  of  the  main  program 
rather  than  in  the  subroutines. 

2.  Elimination  of  the  transmittance  tables  for 
Ho0  vapor,  infrared  0^  and  the  uniformly-mixed 

gases  . 

3.  Reading  of  the  spectral  constants  for  the  newly 
added  band  models  for  the  trace  gases. 

4.  Reading  of  the  air  density  profile  for  the  U.S. 
Standard  atmosphere,  and  of  the  ppm  for  the 
calculation  of  the  equivalent  amounts  of  the 
trace  gases. 

5.  Changes  in  the  dimension  statements  to  include 
the  additional  subscripted  variables. 

Changes  in  the  output  format  include: 

1.  Modification  of  the  print  out  of  the  input  data. 

2.  Modification  of  the  output  table  of  computations 
to  include  the  transmittance  for  the  trace  gases. 

It  should  be  stressed,  however,  that  the  code  is  operated 
using  exactly  the  same  four  control  cards  as  in  the  original 


code. 


8.2  Testing  of  Modularized  Version 

The  first  step  in  the  testing  of  the  modularization 
consisted  of  running  identical  calculations  using  the 
original  code  and  the  modularized  code  before  the  replace¬ 
ment  of  the  transmittance  tables  and  before  the  addition 
of  the  trace  gases.  Numerous  cases  were  considered  during 
this  effort.  A  particular  case  in  which  the  spectral  range 
varied  from  2350  to  2450  cm  ^  for  a  path  at  65°  from  a 
height  of  2.5  km  to  a  height  of  8.5  km  and  a  23  km  visual 
range,  is  shown  in  the  Appendix.  This  output  is  identical 
to  the  output  obtained  from  the  original  Lowtran. 

The  second  step  in  the  testing  procedure  consisted 
of  running  calculations  using  the  original  code  and  the 
modularized  version  with  the  transmittance  tables  replaced 
with  the  continuous  functions,  but  before  the  addition 
of  the  trace  gases.  For  this  purpose,  10  frequencies 
were  selected  such  that  different  combinations  of  models 
would  be  effective  in  the  calculation  of  the  total  trans¬ 
mittance.  The  calculations  were  for  a  5  km  path  at  sea 
level  in  a  sub-arctic  winter  atmosphere  with  a  23  km 
visibility.  The  results  are  summarized  in  Table  9.  The 
columns  listed  under  Transmittance  Deviations  represent 
the  differences  between  the  calculations  using  the  tab¬ 
ulated  and  the  continuous  functions.  Note  that  the  average 
total  transmittance  deviation  is  0.0034,  which  is  below 
the  standard  deviation  obtained  in  the  curve  fitting  of 


TRANSMITTANCE  DEVIATIONS 


WAVENUMBER 
(  cm'  1 ) 


H20 

INFRARED 

UNIFORMLY- 

VAPOR 

MIXED  GASES 

TOTAL 

TRANSMITTANCE 


-J 

455 

0 . 0022 

0.0000 

555 

0.0035 

0.0000 

655 

0.0041 

0.0003 

755 

0.0007 

0.0003 

955 

0.0057 

0.0002 

1155 

0.0026 

0.0003 

1355 

0.0044 

0 . 0000 

1855 

0.0007 

0.0001 

2455 

0.0015 

0.0000 

3155 

0.0037 

0.0001 

0.0000 

0.0021 

0.0026 

0 . 0018 

0.0000 

0 . 0000 

0 .0047 

0.0038 

0.0050 

0.0096 

0.0050 

0 . 0058 

0.0013 

0 . 0006 

0.0034 

0.0007 

0.0054 

0.0045 

0.0027 

0.0053 

rable  9.  Transmittance  difference  between  calculations 
using  the  tabulation  of  the  transmittance  func 
tions  and  calculations  using  the  continuous 
function  representation  for  a  5  km  path  at  sea 
level  in  a  sub-arctic  winter  atmosphere. 
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the  functions  to  the  individual  transmittance  tables.  This 
deviation  amounts  to  an  error  of  about  0.7%  in  the  middle 
of  the  curve-o f -growth ,  which  far  exceeds  the  accuracy  of 
Lowtran  (between  10  to  20%).  The  following  are  attractive 
features  of  the  continuous  functions: 

1.  They  inherently  provide  for  continuous  expo¬ 
nential  interpolation  in  transmittance,  which 
is  superior  to  the  linear  interpolation  used 
in  connection  with  the  transmittance  tables. 

2.  They  provide  for  analytical  operations  such  as 
differentiation  and  interpolation  often  needed 
in  radioactive  transfer  problems. 

3.  They  can  be  used  easily  for  curve  fitting  to 
new  transmittance  data  using  computerized 
procedures . 

4.  Their  use  reduces  significantly  the  computer 
storage  requirements  for  the  individual  models. 

5.  They  continuously  provide  for  transmittance 
calculations  for  small  argument  values  where 
0.9999  <  T  <  1,  for  which  range  Lowtran  4 
includes  an  additional  exponential  function. 

It  should  be  pointed  out  that  the  deviations  listed  in 
Table  9,  although  insignificant,  do  not  represent  errors 
solely  attributed  to  the  analytical  functions.  Since  they 
are  smaller  than  the  uncertainties  in  the  original  data 
used  to  develop  the  tabulated  t r a n sm i t t a nc e s ,  they  pri¬ 
marily  represent  differences  in  the  calculational  procedures. 
In  fact,  in  the  region  between  the  tabulations  the  use  of 
the  analytical  functions  are  likely  to  provide  more  ac¬ 
curate  results  than  the  use  of  the  original  method  in  Lowtran. 


The  last  effort  in  the  testing  of  the  modularized 
code  consisted  of  calculations  involving  the  newly  added 
trace  gases.  For  this  purpose  ten  frequencies  were  run 
at  which  the  trace  gas  models  are  effective.  The  same 
frequencies  were  run  with  the  modularized  Lowtran  without 
these  models.  The  results  are  summarized  in  Table  10. 

The  table  is  primarily  intended  to  show  the  absorptive 


effects  of  the  trace  gases. 


Coefficients  of 

Absorber  Spectral  Parameter  C'  (ci.i"’^)  Analytical  Standard 

arameters  Function  Deviation 


Table  12b,  Band  model  parameters  por  NO 
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Table  12d,  Band  model  parameters  for 


8.3  Band  Model  Development 

Two  sets  of  curves  of  growth  data  for  each  major 
absorption  band  for  four  trace  gases  SC^  NO,  NO^.  and  NH^ 
were  generated  by  the  line-by-line  calculation  from  the 
AFGL  trace  gas  parameter  tape.  One  of  them  consists  of 
12-cut  data  for  several  layers  of  atmosphere  and  the  other 
consists  of  65-cut  data  for  the  standard  atmosphere  only. 
Considering  the  wide  range  of  applications,  we  included 
not  only  the  standard  atmospheric  conditions  but  also  one 
condition  each  from  the  tropical  and  subarctic  winter 
climates.  They  are  listed  in  Table  11  together  with  the 
12  chosen  transmittance  values.  The  major  absorption 
bands  for  the  four  trace  gases  are  given  in  Table  12 
together  with  the  corresponding  computed  C’  values. 

Ten  middle  cuts  were  chosen  from  the  12-cut  data 
and  used  in  both  ADSET  and  SIMMIN  for  the  computation  of 
the  band  model  parameters  and  the  standard  transmission 
function.  Depending  on  the  number  of  major  absorption 
bands,  the  total  numbers  of  data  used  differ  but  are  in 
the  range  of  60-210.  The  65-cut  data  was  used  in  ADSET 
for  the  piecewise  interpolation  to  compute  piecewise 
analytical  transmission  functions. 

The  ADSET  computations  were  done  first.  The  ob¬ 
tained  band  model  parameter  values  n,  m,  and  Cj,  and  nine 
sets  of  coefficients  a^,  a^,  and  a^  are  tabulated  in 
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Table  12  for  SC^,  NO,  NO2,  and  NH^  In  this  order.  The  cor¬ 
responding  standard  deviations  are  also  listed  in  these 
tables  . 

We  also  have  generated  standard  atmospheric  condition 
data  for  non-major  bands  of  each  trace  gases.  These  data 
were  used  to  evaluate  non-major  C'(v)  values.  The  computed 
C ’ (v)  values  were  listed  in  Table  5.  As  we  have  discussed, 
these  C ' (V)  values  and  the  band  model  parameters  together 
with  the  first  order  piecewise-analytical  standard  trans¬ 
mission  function  were  implemented  in  the  modularized 
Lowt  r an . 

We  recall  that  the  SIMM1N  computation  is  a  recursive 
one  and  we  need  a  set  of  initial  guesses  of  the  parameter 
values  to  start  the  computation.  For  the  band  model 
parameters  n,  m,  and  C j ,  we  used  the  values  computed  by 
ADSET.  For  and  a ,,  ,  the  respective  averages  of  the  first 
order  piecewise  interpolation  results  of  ADSET  were  used. 
Finally,  a^  was  set  to  be  zero.  We  note  that  our  initial 
guesses  are  fairly  accurate,  since  these  values  were 
optimal  or  optimal  in  average  for  ADSET  computation.  A 
small  number  £  which  was  used  for  the  check  of  convergence 
was  chosen  to  be  10  ^ .  Since  the  parameter  values  -re 
expected  to  be  in  the  range  -10~10,  £  *  10  ^  gives  the 
limit  of  numerical  accuracy  of  numbers  in  the  computer. 

The  SIMMIN  results  are  also  listed  in  Table  12. 
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Typical  curve-fits  by  piecewise  analytical  standard 
transmission  functions  to  actual  data  are  shown  in  Fig. 

10  for  SO2  at  500  wavenumber.  The  corresponding  analytical 
standard  transmission  function  are  also  compared  to  the 
data  in  Fig.  10.  In  all  of  the  three  graphs  in  this 
Figure,  the  65-cut  data  were  also  plotted  to  show  the 
fitness  of  the  standard  curves. 

The  computation  was  repeated  using  two  smaller 
data  sets  with  6  and  4  cuts  only.  The  chosen  cuts  were 
(0.95,  0.9,  0.8,  0.6,  0.4,  and  0.1)  for  6  cut  data  and 
(0.95,  0.9,  0.6,  and  0.2)  for  4  cut  data.  The  derived 
band  model  parameter  values  were  similar  to  those  in 
Table  12  and,  hence,  were  not  repeated  here.  Instead, 
the  corresponding  standard  deviations  were  listed  and 
compared  with  the  10  cut  cases  in  Table  13. 
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IX  Discussion  and  Conclusions 

9.1  Introduction 

The  modularized  version  presented  here  is  funda¬ 
mentally  the  same  Lowtran  code  except  for  the  separation 
of  its  computation  structure  into  separate  modules  or 
subroutines.  Although  it  is  based  on  the  4th  version,  it 
can  be  adapted  with  little  modification  to  any  future 
versions,  such  as  the  5th  version  now  in  progress.  In 
fact,  this  latter  version  already  has  been  structured  by 
AFGL  such  that  the  emi s s ion / r ad ianc e  loop  is  in  a  subroutine. 
The  modularized  code  presented  here  breaks  down  that  loop 
into  a  frequency  selection  subroutine,  an  equivalent 
absorber  amount  subroutine  and  separate  subroutines  for 
each  one  of  the  attenuation  codes.  The  use  of  modules  in 
a  complex  code  such  as  Lowtran  has  numerous  advantages, 
among  which  the  amenability  for  updating  by  individual 
users  to  suit  their  specific  needs  is  at  the  top  of  the 
list.  In  the  ever  changing  field  of  modeling  it  is  highly 
desirable  to  be  able  to  easily  modify  the  code  for  changes 
in  the  spectral  coverage,  the  spectral  resolution,  the 
absorber  concentrations  in  abnormal  environments,  the 
original  transmission  data  used  in  the  development  and  in 
the  models  used  for  the  individual  attenuators.  The  mod¬ 
ularized  version  presented  here,  although  is  not  the  final 
answer  to  all  conceivable  needs,  it  is  a  first  basic  step 
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in  that  direction.  Practicing  this  predicament,  the 
authors  added  transmission  models  for  the  trace  gases  to 
the  code. 

9.2  Changes  and  Recommendations 

The  following  are  the  basic  changes  introduced  in 
The  Modularized  Lowtran; 

1.  The  original  main  program  was  separated  into  a 
central  program  and  subroutines  for  the  absorber 
amount  and  the  individual  attenuation  models. 

2.  In  the  interest  of  efficiency  and  clarity,  a 
new  subroutine  FGQSL  was  added  for  the  selection 
of  the  attenuation  model  effective  at  the  given 
frequency . 

3.  The  subroutine  HNO^  was  re-structured  to  the 
form  of  the  other  previously  incorporated  sub¬ 
routines  in  Lowtran. 

4.  Continuous  analytical  models  were  provided  to 
replace  the  transmittance  curves  for  H20  vapor, 

and  the  uniformly-mixed  gases. 

5.  New  subroutines  for  the  trace  gases  SO^,  NO,  N02 
and  NH^  were  added. 

A  copy  of  the  modularized  version  is  found  in  the  Appendix. 

Some  recommendations  may  be  made  at  this  time 
concerning  future  modifications  of  Lowtran.  They  are  as 
follows : 

1.  AFGL  should  be  informed  of  the  modularization 
presented  here  as  well  as  of  the  addition  of 
the  trace  gases  (S02>  NO,  N02  and  NH^)  so  that 
they  may  modify  their  master  copy  accordingly. 

2.  As  soon  as  they  are  available,  vertical  profiles 
for  the  concentrations  of  the  trace  gases 
should  be  added. 

3.  The  uniformly-mixed  gases  (C02,  N20,  CO,  CH^ , 
etc.)  should  be  modeled  and  be  included  as 
separate  subroutines. 
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The  resolution  of  all  resonant  absorption  models 
in  the  IR  should  be  increased  to  about  10  cm“l, 
which  will  also  allow  for  model  redevelopments 
with  more  recent  and  more  accurate  transmission 
data . 

All  model  developments  should  adopt  computerized 
numerical  methods  rather  than  the  inaccurate 
manual  graphical  techniques  used  in  the  past. 

The  transmittance  calculations  should  include 
the  calculation  of  the  standard  deviation 
expected  from  known  uncertainties  in  the  input 
meteorological  variables1  . 

The  slant-path  calculations  should  include 
corrections  for  the  Lorentzian-Doppler 
broadening  above  the  10  mb-level. 

Continuous  functions  should  replace  the  tabu¬ 
lated  transmittance  functions,  together  with 
their  awkward  interpolation  procedure. 


9.3  Model  Development 

The  values  of  band  model  parameters  n  and  m  and 
spectral  parameters  Cj  obtained  by  ADSET  and  SIMMIN  agreed 
very  well.  Furthermore,  as  it  was  shown  in  Table  13,  the 
standard  deviations  corresponding  to  different  cases 
followed  a  same  pattern  for  the  ADSET  and  SIMMIN  results. 
This  consistency  proves  the  validity  of  both  methods. 

In  general,  the  SIMMIN  and  ADSET  computations  re¬ 
sulted  in  similar  standard  deviations.  It  was  expected 
that  the  ADSET  computation  should  result  in  lower  standard 
deviations  since  it  contained  more  parameters  to  adjust. 
However,  for  a  half  of  the  cases,  the  SIMMIN  code  produced 
lower  standard  deviations.  This  is  due  to  the  large 
computational  error  for  the  ADSET  computations  in  solving 
the  normal  equation  AX  =  B.  When  the  condition  number  of 
the  coefficient  matrix  A  becomes  large  (i.e.,  A  becomes 
close  to  be  singular),  the  computational  error  becomes  so 
large  that  it  can  exceed  the  directly  minimized  error  of 
the  SIMMIN  computation. 

We  note  that  this  reversal  occurred  for  all  four 
cut  data  cases.  This  suggests  that  the  advantage  for  ADSET 
of  having  more  parameters  to  be  adjusted  is  not  significant 
for  these  cases.  Hence,  we  recommend  the  use  of  SIMMIN 
if  the  available  data  contains  less  than  five  or  six  cuts. 


A  comparison  of  the  standard  deviations  for  two 
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piecewise  interpolation  results  in  the  ADSET  computation 
showed  no  significant  difference.  Furthermore,  the  results 
with  the  second  method  using  quadratic  form  of  x  on  the 
exponent  of  the  double  exponential  function  were  ’bumpy' 
for  some  cases.  Since  the  nature  of  the  transmittance 
does  not  predict  this  behavior,  we  conclude  that  the  first 
method  using  linear  function  of  x  is  accurate  enough  to  be 
used  in  the  actual  application. 

The  standard  deviations  were  much  higher  for  NC^ 
cases  than  the  cases  for  the  rest  of  absorbers.  By  in¬ 
specting  each  curve  of  growth  in  detail,  it  was  found  that 
this  was  mainly  due  to  the  difference  in  the  steepness  of 
the  curves  of  growth  for  three  absorption  bands.  This 
difference  cannot  be  compensated  by  C|  values  since  they 
only  shift  the  curves  of  growth  linearly.  In  fact,  within 
the  current  band  model  structure,  it  is  impossible  to 
compensate  this  difference.  Hence,  it  may  be  necessary  to 
modify  some  of  the  basic  assumptions  regarding  the  band 
model  structure,  if  lower  standard  deviations  are  required. 

As  a  side-effect  of  this  discrepancy  in  the  tangent 
of  curves  of  growth,  the  SIMMIN  computation  took  far  more 
time  for  NO^  cases  than  the  rest.  Most  of  the  computations 
of  ADSET  were  completed  by  26-36  CPU  seconds.  The  fluc¬ 
tuations  in  the  computation  time  were  very  small.  On  the 
other  hand,  the  SIMMIN  computation  time  varied  from  14 
seconds  to  270  seconds.  NO2  cases  consumed  about  200-270 


seconds,  which  were  about  four  times  as  much  as  that  for 
the  other  cases.  This  is  because  the  minimizing  point  in 
the  parameter  space  is  not  well  defined  for  NO2  cases.  In 
other  words,  the  error  surface  in  the  parameter  space  has 
a  very  shallow  bottom  so  that  the  updating  step  cannot 
produce  large  enough  changes  in  the  parameter  guesses  in 
order  to  have  a  rapid  convergence. 

Thus,  it  was  found  that  the  accuracy  of  the  computed 
results  and  the  time  of  execution  depend  heavily  on  the 
actual  data.  Hence,  it  is  very  important  to  give  enough 
consideration  for  the  data  structure.  This  will  be  dis¬ 
cussed  in  the  next  section. 
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9.4  Data  Structure 

As  it  was  expressed  earlier,  we  assumed  that  the 
number  of  layers  (=  the  number  of  data  points)  in  each  cut 
is  the  same  for  an  absorption  band.  This  was  done  for  the 
sake  of  easier  coding  in  data  handling.  However,  this 
assumption  need  not  be  valid.  Especially  in  weaker  absorp¬ 
tion  bands,  it  is  required  to  use  very  large  range  values 
to  have  high  enough  equivalent  absorber  amounts  in  order 
to  realize  lower  t r ansmi t tanc es .  In  some  cases  the  range 
becomes  enormous  (in  the  order  of  the  radius  of  the  earth) 
so  that  the  corresponding  data  no  longer  possess  physical 
significance.  The  ADSET  code  has  a  criterion  that  if  the 
logarithm  of  the  equivalent  absorber,  log  W,  exceeds  a 
certain  critical  value,  then  the  corresponding  data  will 
be  set  aside  and  will  not  be  used  in  the  later  computation. 
The  critical  value  was  set  to  be  2  for  the  actual  computa¬ 
tion,  which  corresponds  approximately  to  a  vertical  path 
through  the  atmosphere. 

In  connection  with  this,  if  data  are  not  available 
at  some  layers,  then  the  data  values  are  set  at  0  to  flag 
the  nonavailability  of  data.  ADSET  can  also  detect  this 
and  will  ignore  the  data. 

A  caution  must  be  executed  in  choosing  combinations 
of  pressures  and  temperatures,  i.e.,  atmospheric  conditions. 
If  a  data  set  contains  either  the  standard  pressure  or  the 


standard  temperature  or  both  only,  then  both  ADSET  and 
SIMMIN  fail  because  of  the  fact  that  the  coefficient  of  n 
or  m  or  both  in  Eq.(12)  becomes  zero,  since 

P 

log  (~)  -  log  1  -  0,  (61) 

o 

T 

log  =  log  1  =  0.  (62) 

o 

For  this  case,  the  coefficient  matrix  A  of  the  normal 
equation  in  ADSET  becomes  singular  and  the  gradient  cor¬ 
responding  to  n  or  m  or  both  in  SIMMIN  becomes  zero  all  the 
time.  Hence,  the  normal  equation  cannot  be  solved  in  ADSET 
and  the  initial  guess  of  n  or  m  or  both  cannot  be  changed 
in  SIMMIN. 

Another  consideration  which  should  be  pointed  out 
is  to  include  different  climate  conditions.  The  standard 
climate  condition  for  several  layers  of  atmosphere  contains 
sequence  of  pressures  and  temperatures  both  of  which  are 
monotone  decreasing.  Therefore,  if  only  these  conditions 
are  used,  then  it  is  very  difficult  to  distinguish  the 
cause  of  changes  in  the  transmittance  due  to  the  changes 
in  pressure  and  in  temperature.  This  leads  to  the  shallow 
bottom  of  the  error  surface  and  hence,  large  computational 
error  results  in  ADSET  caused  by  the  large  condition 
number  of  the  coefficient  matrix  and  slow  convergence  in 
SIMMIN  due  to  the  small  gradient.  In  the  actual  computation 
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we  Included  not  only  the  standard  climate  conditions  but 
also  one  condition  each  from  the  tropical  and  subarctic 
winter  climates  in  consideration  of  wide  applicability  of 
the  results.  Numerically  speaking,  this  also  resulted  in 
making  the  regression  problem  well-posed  by  breaking  the 
monotonousness  of  the  pressure  and  temperature  combinations 
of  the  standard  conditions.  In  fact,  several  computations 
were  done  for  AOSET  and  SIMMIN  with  standard  condition  data 
only.  SIMMIN  took  10-45  minutes  of  CPU  time  to  converge 
if  it  were  convergent  and  ADSET  resulted  in  a  set  of  absurd 
values  for  n  and  m.  Thus,  the  importance  of  the  numerical 
consideration,  which  is  ignored  in  many  cases,  is  clearly 
indicated.  The  proper  care  should  be  taken  when  selecting 
controllable  data  values. 


A 
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ADSET  V  1-3 

EVALUATION  OF  ABSORBER  PARAMETERS  AND  STANDARD  EMPIRICAL  AND 
PIECEWISE -ANALYTICAL  TRANSMISSION  FUNCTIONS 

THIS  CODE  USES  THE  SUBROUTINE  SIMQ  IN  SSP  LIBRARY 

THIS  CODE  CONSISTS  OF 

MAIN:  COMPUTATIONS  OF  BAND  PARAMETERS  N,M,C' 

AND  EMPIRICAL  STANDARD  TRANSMISSION  FUNCTION 
NMBC :  COMPUTATION  OF  NON-MAJOR  BANDS'  C-VALUES 

INTPL1:  COMPUTATION  OF  THE  PIECEWISE- ANALYTICAL 

STANDARD  TRANSMISSION  FUNCTION  GIVEN  BY 
TAU  =  EXP(-10**(A1+A2*X)  ) 

INTPL2 :  COMPUTATION  OF  THE  PIECEWISE- ANALYTICAL 

STANDARD  TRANSMISSION  FUNCTION  GIVEN  BY 
TAU  =  EXP(-10**(A1+A2*X+A3*X**2)) 

SDTAU :  COMPUTATION  OF  THE  ERROR  STANDARD  DERIVATIONS 

BETWEEN  PIECEWISE  ANALYTICAL  STANDARD 
TRANSMISSION  FUNCTION  AND  THE  ORIGINAL  DATA 
USED  IN  THE  MAIN  PROGRAMME 

DATA  SET-UP 

1 .  1ST  CARD  TITLE  IN  20A4 

(ABSORBER  TYPE, ETC) 

2.  2ND  CARD  FOUR  CONTROL  NUMBERS  IN  415 

MAXRPT  MAXIMUN  NUMBER  OF  REPETITION  OF 
THE  COMPUTATION  IN  MAIN 
INDX ( 1 )  SUBROUTINE  NMBC  IS  CALLED  IF  >  0 
INDXC2)  SUBROUTINE  INTPL1  IS  CALLED  IF  >  0 
INDX ( 3 )  SUBROUTINE  INTPL2  IS  CALLED  IF  >  0 

3.  DATA  SET  FOR  MAIN 

CONSISTS  OF  SEVERAL  SUBSETS  (MAX. 6)  OF  DATA 
CORRESPONDING  TO  INDIVIDUAL  MAJOR  BANDS 
FIRST  CARD  FOR  EACH  SUBSET  IS  A  CONTROL  CARD  WHICH 
CONTAINS  BAND//,  WAVE  NUMBER,  //  OF  CUTS  AND  //  OF  LEVELS 
IN  THIS  ORDER  BY  THE  FORMAT  ( 1 5 , F 1 0 . 3 , 21 5 ) 

REFER  TO  THE  READ(5,105)  AND  FORMAT  105  FOR  THE 
CONTENTS  OF  EACH  CARD. 

END  OF  DATA  IS  DEFINED  BY  A  BLANK  CONTROL  CARD 

4.  DATA  SET  FOR  NMBC 

INDIVIDUAL  DATA  FORMAT  -  SAME  AS  FOR  MAIN 
END  OF  EACH  DATA  SET  FOR  A  BAND  IS  MARKED  BY 
A  BLANK  CARD 

END  OF  ALL  DATA  IS  MARKED  BY  -1  (12)  IN  ADDITION  TO 
A  BLANK  CARD 

IF  NO  DATA  BUT  A  BLANK  CARD  IS  SUPPLIED,  THEN 
THIS  SUBROUTINE  IS  SKIPPED. 

5.  DATA  SET  FOR  INTPL1 

FORMATION  OF  THE  DATA  IS  THE  SAME  AS  THAT  FOR  MAIN 
DATA  (IF  SUPPLIED)  WILL  BE  USED  FOR  S.D.  COMPUTATION 
ONLY. 

IF  NO  DATA  BUT  A  BLANK  CARD  IS  SUPPLIED,  THEN  THE 
S.D.  COMPUTATION  IS  SKIPPED. 

6.  DATA  SET  FOR  INTPL2 
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FORMATION  OF  THE  DATA  IS  THE  SAME  AS  THAT  FOR  MAIN 
IF  NO  DATA  BUT  A  BLANK  CARD  IS  SUPPLIED,  THEN 
THIS  SUBROUTINE  IS  SKIPPED. 

NOTE:  DATA  SET  MUST  HAVE  A  CUT  STRUCTURE  SUCH  THAT  EQUAL 

TRANSMITTANCE  DATA  ARE  GROUPED  TOGETHER  AND  THESE  GROUPS 
ARE  QUEUED  IN  THE  DECENDING  ORDER  IN  TAU.  THE  QUEUING  OF 
THE  LEVELS  WITHIN  EVERY  GROUP  MUST  BE  THE  SAME. 

DIMENSION  V(19) ,A(19, 19) ,X(361 ) ,B(19) ,RI(6, 12, 10) 

DIMENSION  POO)  ,  WWW  (12, 1 0 )  ,  STANDV ( 1 2 ) , TSD( 6 ) , NDATAC 6 ) , INDX ( 3 ) 
COMMON  /PARM  1  /  TSTD ( 1 2 ) , PW( 1 2 ) , WN ( 6 ) , CSTD (6 ) , NCUT, NC , NAME (20 ) , 

*  AN, AM, CF, ICON ST (6) ,NEL 

COMMON  /PARM2/  PRES (6 , 1 2 , 1 0 ) , TEMP (6 , 1 2 , 1 0 ) , UGAS (6 , 1 2 , 1 0 ) , 

*  TAU(6,12),NTC(6), NLV (6 ) 

CF  =  1 .0 

LOOPCT=  1 
WCRIT=2. 

M  =  0 

READ(5 ,100)  (NAME(I) ,1=1 ,20) 

100  FORMAT (20A4 ) 

READ(5, 101 )  MAXRPT, (INDX(I) ,1=1 ,3) 

101  FORMAT (415 ) 

COMPUTATION  OF  ABSORBER  PARAMETERS  N,  M  &  C-VALUES  IS  REPEATED 
MAXRPT  TIMES,  WHERE  1  <  MAXRPT  <  10  IS  READ  IN  BY  15  FORMAT 
(SUGGESTED  VALUE  IS  5) 

DATA  READ-IN  ROUTINE 

1000  CONTINUE 

READ(5 ,102)  IC , W , JM , KM 

102  F0RMAT(I5,F10. 3,215) 

IF(IC.LE.O)  GO  TO  2000 
IF(M.GT.O)  GO  TO  10 

CALL  DATE  ( MONTH , IDAY , IYEAR ) 

WRITE (6, 1 1 1 ) MONTH, IDAY, IYEAR 
111  FORMAT( 1H1 ,T60, 14, '  /  ',12,'  /  ’,12,//) 

WRITE ( 6 , 200 )  (NAME(I),I=1 ,20) 

200  FORMAT ( 1 H  ,T25,20A4) 

GO  TO  11 

10  CONTINUE 
WRITE(6 ,201  ) 

201  FORMAT ( 1 H 1 ) 

11  CONTINUE 
M=M+1 

WN ( M) =W 
NTC ( M) = JM 
NLV ( M) =KM 

WRITE (6, 202)  M,WN(M) ,NTC(M) ,NLV(M) 

202  FORMAT ( 1  HO , T 1 5 , ' ***  BAND', 13,'  ( WAVE  NUMBER  =', F 1 0 . 3 ,' ) 

*  ///,T20, 'TOTAL  #  OF  CUTS  =', 13 ,//, T20 , 'TOTAL  #  OF  LEVELS  =’,I3, 

*  ///) 

WRITEC6 ,203) 

203  FORMAT ( 1H  ,T5,'(  DATA  FORMAT  )',//, T9 , 'GAS# ', T 15 ,' WAVE  #’,T24, 

*  'PRESSURE' ,T36, 'TEMP. ’ ,T46, 'PPM' ,T58, 'RANGE' ,T70, 'UGAS' ,T78, 

•  'TRANSM.',/) 
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DO  12  J=1 , JM 
WRITE ( 6 , 204  )  J 

204  FORMAT ( 1H0 , T5 , ' <  CUT', 13,'  >’,/) 

T =0 . 

IT  =  0 

DO  13  K  =  1 , KM 

READ(5, 103)  KGAS , FREQ , RPRES , RTEMP , PPM , RANGE , RUGAS , TX 
103  F0RMAT(I2,F10.3,E11 .4,F9.3,E11.4fE13.6,E11.4,F7.4) 

RUGAS=RUGAS/CF 

WRITE (6, 205)  KGAS , FREQ , RPRES , RTEMP , PPM , RANGE , RUGAS , TX 

205  F0RMAT(T10,I2,F10.3,E11.4,F9.3,E11 .4, El  3. 6, Ell .4.F8.4) 

PRES,  TEMP  &  UGAS  ARE  CONVERTED  TO  THE  LOG  OF  THE  NORMALIZED 
VALUES.  IF  RPRES=0  (INDICATES  NO  DATA), THEN  UGAS(M,J,K)  IS  SET 
AT  AN  IMPOSSIBLE  VALUE  ,  ALSO  RI(M,J,K)  IS  SET  TO  ZERO. 

IF ( R PRES . GT . 0 . )  GO  TO  14 
PRES ( M , J , K )  =  0  . 

TEMP (M , J , K) =0 . 

UGAS ( M ,  J  ,  K )  =  1  0  . 

RI(M, J,K)sO. 

GO  TO  15 

14  CONTINUE 

PRES(M, J,K)=ALOG10(RPRES/1013. ) 

TEMP ( M , J , K ) = ALOG 10(273- 15/RTEMP) 

UGAS(M,J,K)=ALOG 10 (RUGAS) 

15  CONTINUE 

T=T+TX 
IT=IT+1 
RI ( M , J , K ) = 1 .0 
13  CONTINUE 

TAU ( M , J ) =T/ FLOAT ( IT ) 

12  CONTINUE 
GO  TO  1000 

END  OF  DATA  INPUT 

CONSTANTS  USED  IN  LATER  COMPUTATION  ARE  INITIALIZED 
FROM  2000  TO  3000. 

NCUT  =  MAXIMUM  //  OF  CUTS  USED  IN  COMPUTATION 
NDIM  =  DIMENSION  OF  THE  COEFFICIENT  MATRIX 

2000  CONTINUE 

IF(M.GT.O)  GO  TO  20 
WRITE (6 , 206  ) 

206  FORMAT(1HO,//,T10, '$$$  NO  INPUT  DATA  $$$') 

STOP 

20  CONTINUE 
NC=M 

CSTD(1 )=0 . 

NCUT  =  NTC ( 1  ) 

NPTS=NTC ( 1 ) *NLV( 1 ) 

IF ( NC . LE . 1 )  GO  TO  21 
DO  22  1=2, NC 
NCUT=MAXO (NCUT , NTC ( I ) ) 
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NPTS=NPTS+NTC(I)*NLV(I) 

22  CONTINUE 
21  CONTINUE 

FNC=FLOAT (NC) 

RIT=FLOAT (NPTS) 

DO  23  J  =  1  ,  NC UT 
TC  =0  . 

DO  24  M= 1 , NC 
TC=TC+TAU(M, J) 

24  CONTINUE 

TSTD ( J ) =TC/FNC 

23  CONTINUE 
NDIM=NC+1 +NCUT 

t*«**«««««******«****»****«*«*«****«*»*««****»*»««***«*******««**«**ft** 

COMPUTATION  OF  THE  ABSORBER  PARAMETERS. 

THIS  LOOP  WILL  BE  REPEATED  MAXRPT  TIMES. 

FORMATION  OF  THE  NORMAL  EQUATION  AX  =  B  ,  A  IS  SYMMETRIC 

3000  CONTINUE 

DO  30  1=1,19 
B  ( I )  =0  . 

DO  31  J= 1 , 1 9 
A ( I , J ) =0 . 

31  CONTINUE 
30  CONTINUE 

DO  1  M=1 , NC 
JM  =NTC ( M ) 

KM=NLV(M) 

DO  2  J=1 , JM 
DO  3  K  =  1 , KM 

IF(RI(M, J,K) .LT.0.5. )  GO  TO  3 

DO  4  IC  =  1  ,  19 

V(IC)=0. 

4  CONTINUE 

V  (NDIM-M )  =  1 . 

V(NDIM)=PRES(M,J,K) 

V (NDIM-1 ) =TEMP ( M , J , K  ) 

VV=-UGAS ( M , J , K ) 

V (NCUT+1 -J  )  =  1 . 

DO  5  11=1 , NDIM 
DO  6  IJ=1 , NDIM 

A(II,IJ)=V(II)*V(IJ)*RI(M,J,K)  +  A ( I I , I J ) 

6  CONTINUE 

B(II)=V(II) *VV*RI ( M , J , K)  +  B  ( 1 1 ) 

5  CONTINUE 
3  CONTINUE 
2  CONTINUE 
1  CONTINUE 

IF  J-TH  ROW  OF  "A"  IS  ZERO,  A(J,J)  IS  CHANGED  TO  -1 
WHICH  IS  DONE  IN  ORDER  TO  MAKE  "A”  NON-SINGULAR 
THIS  HAPPENS  WHEN  ALL  OF  THE  DATA  FOR  BAND  J-1  FAIL  TO  SATISFY 
THE  CRITERION  W  <  WCRIT.  THE  BAND  J-1  WILL  BE  IGNORED  IF  THIS 
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HAPPENS,  AND  THE  C-VALUE  FOR  BAND  J-1  WILL  BE  COMPUTED 
C  SEPARETELY. 

C 

ICONST ( 1 )=1 
IF ( NC . EQ . 1  )  GO  TO  MO 
DO  Ml  M=2 , NC 
ICONST (M)  =  1 
I =NDIM-M 

IF ( A( I , I ) . NE . 0 . )  GO  TO  Ml 
A  ( I ,  I )  =  - 1 .0 
ICONST ( M) =0 
Ml  CONTINUE 
MO  CONTINUE 
C 
C 

NCOL=0 

DO  M2  J=1 , NDIM 
DO  M3  1=1 , NDIM 
NCOL=NCOL+ 1 
X(NCOL)=A( I , J ) 

M3  CONTINUE 
M2  CONTINUE 

PRINTING  OF  THE  HEADING  FOR  EACH  TRIAL  AND  THE  NORMAL  EQUATION 

IF ( LOOPCT . GT . 1 )  GO  TO  50 
WRITE ( 6 , 207 )  MAXR PT , LOOPCT 

20  7  F0RMAT(1H1  ,  T20  ,  '***  ABSORBER  PARAMETER  COMPUTATION  *••',///, 

*  T 15 , 'NOTE :  THE  COMPUTATION  WILL  BE  REPEATED  MAXRPT  s', 12, 

*  '  TIMES. ',////, T10, 'TRIAL  # ' , II , 5X , » ( ALL  DATA  WERE  USED)’) 

GO  TO  51 

50  CONTINUE 
WRITE (6 , 208 )  LOOPCT 

208  F0RMATC1H1 , T 1 0 , 'TRIAL  1 1 , 5X ,' (PARTIAL  DATA  WERE  USED  WITH’, 

*  '  CUT-OFF  CRITERION  :  W  <  2  )’) 

51  CONTINUE 

WRITE(6 , 209 )  NDIM, NDIM 

209  FORMAT ( // , 1  HO , ' <  NORMAL  EQUATION  :  AX  =  B  >',//, T 1 0 ,', WHERE  THE', 

*  '  COEFFICIENT  MATRIX  AC, 13,'  , ',13,'  )  AND  THE  CONSTANT  VECTOR', 

*  ’  B  ARE ’ ,//) 

IFCNDIM.LE. 17)  GO  TO  52 
WRITE (6 ,210)  NDIM 

210  FORMAT ( 1 H  ,'***  WARNING  :  DIMENSION  OF  THE  MATRIX  IS  TOO  LARGE', 

*  '  (',13,'  )  TO  BE  PRINTED  IN  A  MATRIX  FORM  •*»',/) 

52  CONTINUE 
DO  53  1=1 , NDIM 

WRITE (6 ,21 1 )  ( A ( I , J ) , J  =  1 , NDIM)  ,B(I ) 

211  FORMAT ( 1H  ,  1 8F7 - 3 ) 

53  CONTINUE 

*****  MATRIX  INVERSION  SUBROUTINE  SIMQ  IN  SSP  IS  CALLED  **«*«*»*•«»**• 
CALL  SIMQ(X,B,NDIM,KS) 

PRINTING  OF  THE  SOLUTION  FOR  THE  NORMAL  EQUATION 
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IF(KS.EQ.I)  WRITE (6, 212) 

212  F0RMATOH0.T10, 'WARNING:  THE  COEFFICIENT  MATRIX  IS  SINGULAR.') 
AN=B(NDIM) 

AM=B ( NDIM-1 ) 

IF ( NC . LE . 1 )  GO  TO  54 
DO  55  M=2 , NC 
CSTD(M)=B(NDIM-M) 

55  CONTINUE 
54  CONTINUE 

DO  56  J=1 , NCUT 
PW(J)=-B(NCUT+1-J) 

56  CONTINUE 

WRITE (6, 21 3)  AN, AM, (CSTD(M) ,M=1 ,NC) 

213  FORMATC// ,  1H0, '  <  RESULTS  >',///, T7 N T1 7 M T27 C 1 T37 C2 ’ , 

*  T47, 'C3' ,T57, 'C4' ,T67, 'C5' ,T77, 'C6' ,//,2F10.5,6F10.3) 

WRITE (6,214)  (PW(I) , 1=1 , NCUT) 

214  FORMAT (/ , 1H0 , T7 , ' X *  1 ’ , T 1 7 , ' X*2 ' , T27 , 'X*3' ,T37, 'X«4' ,T47, ' X*5 ' , 

»  T57,  'X*6' ,T67, 'X*7 ' ,T77, 'X»8 ' ,T87, 'X*9' ,T97, ’X*10’ ,T107, 'X«1 1 ' , 

«  T 1 1 7 , ' X *  1 2 ' ,//, 12F10.3) 

NEL=NPTS-INT(RIT) 

WRITE (6 ,215)  NEL 

215  FORMAT(// , 1H0,T4, '#  OF  ELIMINATED  POINTS  =',I5) 

CHECKING  OF  THE  CRITERION  (  W  <  WCRIT  )  AND  THE  COMPUTATION 
OF  C-VALUES  FOR  THE  IGNORED  BANDS. 

RI(M,J,K)  =  0  IF  W  IS  GREATER  THAN  OR  EQUAL  TO  WCRIT 
R I ( M , J , K )  =  1  IF  W  IS  LESS  THAN  WCRIT 

RIT=0 . 

DO  60  M= 1 , NC 
JM  =NTC ( M) 

KM=NLV ( M) 

C AVG=0 . 0 
DO  61  J=1 , JM 
DO  62  K  =  1 , KM 

W=AN*PRES(M,J,K)+AM*TEMP(M,J,K)+UGAS(M,J,K) 

IFCW.GE .WCRIT)  RI(M,J,K)=0. 

RIT=RIT+RI(M,J,K) 

CAVG=CAVG+(PW(J)-W) 

62  CONTINUE 
61  CONTINUE 

IF( ICONST(M) . EQ. 1 )  GO  TO  60 
CSTD(M)=CAVG/FLOAT(JM*KM) 

WRITE (6 , 2 1 6 )  M , M , M , CSTD( M) 

216  FORMAT(/// ,  1H  ,T7,'»*  WARNING  *»’,T25,'NO  DATA  FOR  BAND’, 12, 

*  '  SATISFIES  THE  CRITERION  (  W  <  2  ).',//, T25 THE  C',11, 

*  '  VALUE  IS  SEPARATELY  COMPUTED  BY  AVERAGING .',//, T30 ,' C ', 1 1 , 

*  '  =  '  ,  F 1 0 . 3  ) 

60  CONTINUE 

COMPUTATIONS  OF  STANDARD  DEVIATIONS  IN  X 

NGDATA=0 
GTSD=0 . 

ICST=NC 
DO  70  M  =  1 , NC 
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JM=NTC ( M) 

KM=NLV(M) 

NDATA( M)  =  0 
TSD(M) =0 . 

WRITE ( 6 , 20 1 ) 

WRITE (6 ,202)  M , WN (M) , NTC ( M) , NLV(M) 

WRITE (6 ,217)  AN , AM , M , CSTD(M) 

217  FORMAT(1H  ,T10,'N  = ' , F 1 0 . 5 , // , T 1 0 , ' M  = • , F 1 0 . 5 , // , T 1 0 , ' C ' , 1 1 , 

*  '  =  '  ,  F 1 0 . 5 ) 

WRITEC6 ,218) 

218  FORMAT(//, 1H0,T7, ’RECOMPUTED  X-VALUES  AND  STANDARD  DEVIATIONS', 

*  '  IN  X-VALUES' ,/, 1H0,T2, ’CUT' , Til , 'TAU' ,T20, ’X«' ,T30, 'XI ’ ,T39, 

*  'X2' ,T48, 'X3' ,T57, 'X4' ,T66, 'X5' ,T75, 'X6' ,T84, 'X7' ,T93, 'X8' , 

*  T102, 'X9' , T 1 1 1 , ' X 1 0 • , T 1 2 1 , 'CUTWISE-SD' ,/) 

COMPUTATION  OF  THE  CUTWISE  STANDARD  DEVIATIONS  IN  X 

DO  71  J  =  1  ,  JM 
DN  =  0  . 

WW=0. 

DO  72  K= 1 , KM 

P(K)=CSTD(M)+AN*PRES(M,J,K)+AM#TEMP(M,J,K)+UGAS(M,J,K) 

WWW( J ,K)= (PW(J)-P(K) )*»2*RI(M, J,K) 

WW=WW+WWW ( J , K ) 

DN=DN+RI ( M , J , K ) 

72  CONTINUE 
WW=SQRT ( WW/DN ) 

NDATA(M)=NDATA(M)+IFIX(DN) 

WRITE (6 ,219)  J,TAU(M,J),PW(J) ,(P(K)  ,  K  =  1 ,KM) 

219  FORMAT ( 1 H  , 15 ,F9 . 3 ,F9 . 4 , IX , 10F9 . 4) 

WRITE (6 , 220 )  WW 

220  FORMAT(1H+,T121 , F 1 0 . 5 ) 

71  CONTINUE 

COMPUTATION  OF  THE  LEVELWISE  STANDARD  DEVIATIONS  IN  X 

DO  73  K  =  1 , KM 
WW  =  0. 

DN  =0  . 

DO  74  J=1 , JM 
WW=WW+WWW(J ,K) 

DN  =DN+R I ( M , J , K ) 

74  CONTINUE 

TSD( M) =TSD(M)+WW 
ST ANDV(K)=SQRT( WW/DN) 

73  CONTINUE 

WRITE ( 6 , 22 1 )  (STAND V(K) ,K=1 ,KM) 

221  F0RMAT(1H0,T4, 'LEVELWISE-SD  : ’ , T26 , 1 0F9 . 5 ) 
NGDATA=NGDATA+NDATA(M)*ICONST(M) 

GTSD=GTSD+TSD( M) *FLOAT( I CONST ( M) ) 

ICST=ICST-ICONST(M) 

TSD(M)=SQRT(TSD(M) /FLOAT (NDATA(M) ) ) 

WRITE ( 6 , 222 )  TSD(M) 

222  FORMAT(// , 1H0,T4, 'TOTAL  STANDARD  DEVIATION  FOR  THIS  BAND  :', 

»  F15.6) 

70  CONTINUE 
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PRINTOUT  OF  THE  SUMMARY. 

ALL  VITAL  INFORMATIONS  ARE  PRINTED  OUT  HERE. 


GTSD=SQRT (GTSD/FLOAT (NGDATA ) ) 

WRITE (6 , 223 )  LOOPCT , AN , AM 

223  F0RMATOH1  ,  T15  ,  '  ***  SUMMARY  OF  THE  ABSORBER  PARAMETER’, 

*  '  COMPUTATION  FOR  TRIAL  #', 12,'  *«*',///, T20 , 

*  'PRESSURE  EXPONENT  N  = ' , F 1 0 . 5 , // , T20 , 

*  'TEMPERATURE  EXPONENT  M  = ’ ,F 1 0 . 5 , /// , T5 , ' CASE  # ' , 3X , 

*  'WAVE  NUMBER' ,5X, 'C-VALUE' ,5X, 'TOTAL  #  OF  DATA',3X, 

*  'CASEWISE  S.D.  IN  P') 

WRITE (6, 224) * (M,WN(M) ,CSTD(M) , NDATA( M) , TSD(M) , M= 1 , NC ) 

224  FORMAT(1HO,T6,I3,6X,F9.2,5X,F8.3, 10X,I3, 12X.F12.6) 

WRITE ( 6 , 225 )  NGDATA , NEL , GTSD 

225  FORMATC//, 1H0.T15, 'GRAND  TOTAL  #  OF  DATA  = ' , 15 , // , T 1 5 , ' #  OF', 

*  '  ELIMINATED  DATA  =', 15, // , T 1 5 , 'GLOBAL  STANDARD  DEVIATION  IN  P', 

*  '  ='  ,  F 1 2 . 6 , //  ) 

IF(ICST.LE.O)  GO  TO  75 
DO  76  M= 1 , NC 

IF ( ICONST (M) . EQ . 1 )  GO  TO  76 
WRITE(6 , 226  )  M 

226  FORMAT ( 1 H  ,T15,'NOTE:  THE  BAND', 13,'  IS  NOT  INCLUDED  IN  THE', 

«  '  FINAL  STANDARD  DEVIATION’) 

76  CONTINUE 
75  CONTINUE 

WRITE (6, 227)  LOOPCT , (TSTD(J ) , PW ( J ) , J= 1 , NCUT ) 

227  FORMAT(///, 1H 0 , T 1 5 , '***  STANDARD  EMPIRICAL  TRANSMISSION', 

*  '  FUNCTION  FOR  TRIAL  #',I2,'  ***’,///,  T20  ,'  TAU  ',  T35  ,’  X*  ’,/ , 

*  (1HO,T17,F7.3,T30,F8.4)) 


IFCRIT.GT.O. )  GO  TO  80 

IF  NO  INPUT  DATA  SATISFIES  THE  CRITERION,  THE  COMPUTATION  IS 
TERMINATED.  THE  MOST  RECENT  RESULTS  WILL  BE  USED  IN  THE  SEQUAL 

WRITE(6 ,228) 

228  F0RMAT(1H1,//,T15, '$$$  NO  INPUT  DATA  SATISFIES  THE  CRITERION  OF' 

*  '  (  W  <  2  )  $$$' ,//,T15, '$$$  THE  COMPUTATION  FOR  THIS  STEP  IS' 

*  '  TERMINATED  $$$') 

GO  TO  4000 

80  CONTINUE 

LOOPCT  =  L00  PCT  + 1 

IF ( LOOPCT . GT . MAXR  PT )  GO  TO  4000 
GO  TO  3000 

4000  CONTINUE 

SUBROUTINE  COMPUTATIONS  FOLLOW 

IF ( I NDX ( 1 ) .LE.O)  GO  TO  90 
C 

CALL  NMBC 
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90  CONTINUE 

IF(INDX(2) . LE . 0 )  GO  TO  91 
C 

CALL  INTPL1 
C 

91  CONTINUE 

IF ( INDX ( 3 ) . LE . 0 )  GO  TO  92 
C 

CALL  INTPL2 

92  CONTINUE 
STOP 
END 

SUBROUTINE  NMBC 

COMPUTATION  OF  C'-VALUES  FOR  NON-MAJOR  BANDS 
DIMENSION  B( 15) ,CS( 15)  ,FS( 15) 

COMMON  /PARM 1 /  TSTD ( 1 2 ) , PW ( 1 2 ) , WN ( 6 ) , CSTD ( 6 ) , NCUT , NC , NAME (20 )  , 

*  AN, AM,CF, ICONST(6) ,NEL 
WRITE(fi,5)  (NAME(I) ,1=1 ,20) 

5  F0RMAT(1H1 ,T15,20A4) 

WRITE (6 , 10) 

10  FORMAT ( 1  HO , T 1 5 , '  ***  CALCULATION  OF  THE  SPECTRAL  PARAMETERS', 

*  'FOR  NON-MAJOR  BANDS  ***'///) 

DF  =  1 . E30 

11  CONTINUE 
NFREQ=0 

12  CONTINUE 
C  =  0. 

1=0 

15  CONTINUE 

READ(5 , 20 )  KGAS,FREQ,P,T,UGAS,TX 
20  FORMAT(I2,F10.3,E1  1  . 4 ,F9 . 3 , 24X, El  1 . 4 ,F7 . 4) 

IF(KGAS.EQ.O)  GO  TO  25 
IF  (KGAS.LT.O)  GO  TO  35 

THE  FOLLOWING  IF -STATEMENT  IS  INSERTED  TO  DETECT 
AND  TO  IGNORE  THE  INVALID  DATA  POINTS. 

IF(UGAS.GE.DF)  GO  TO  15 
1  =  1  +  1 
WX=FREQ 
UGAS=UGAS/CF 

C=C+(PW(I )-AN*ALOG10(P/1013. ) -AM* ALOG 1 0 ( 27 3 . 15/T)-ALOG10 (UGAS) ) 
GO  TO  15 
25  C=C/FLOAT ( I ) 

NFREQ=NFREQ+1 
CS ( NFR  EQ )  =  C 
FS (NFREQ) =WX 
DO  27  M  =  1 , NC 

IF ( ABS ( WX-WN ( M) ) .LE.O. 1 )  CS ( NF REO ) =CSTD ( M) 

27  CONTINUE 

IF(NFREQ. EQ. 10)  GO  TO  30 
GO  TO  12 
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30  CONTINUE 

WRITE (6,31)  ( FS (K ) , K  =  1 , NFREQ) 

31  FORMAT(1HO,2X, 'WAVE  NUMBER 2X , 1  OF  1 1 . 0 ) 

WRITE (6 , 32 )  ( CS (K ) , K  =  1 , NFREQ) 

32  FORMAT ( 1  HO , 5X , ' C  VALUES ' , 2X , 1  OF  1  1  .  3// ) 

GO  TO  11 

35  CONTINUE 

IF (NFREQ.  EQ.  0)  GO  TO  40 
WRITE (6 ,31)  (FS(K) ,K=1 , NFREQ) 

WRITE (6 , 32 )  ( CS (K ) , K  = 1 , NFREQ) 

40  CONTINUE 
RETURN 
END 

SUBROUTINE  INTPL1 

COMPUTATION  OF  THE  STANDARD  PIECEWISE-ANALYTICAL  TRANSMISSION 
FUNCTION 

VERSION  1-1  **  A3 ( I )  =  0  ** 

TAU  =  EXP ( - 1 0 ** (  A 1 ( I ) +A2 ( I )  *X  )) 

DIMENSION  SDCUT ( 1 5 ) , IC UT ( 1 5 ) ,SDTCUT(15) , ITCUTC 1 5 ) 

COMMON  /P ARM  1 /  TSTD ( 1 2 ) , PW ( 1 2 ) , WN ( 6 ) , CSTD ( 6 )  ,  NC UT , NC , NAME ( 20 ) , 

*  AN, AM, CF, IC0NST(6) ,NEL 

COMMON  /PARM3/  A 1 ( 1 1 ) , A2 ( 1 1 ) , A3 ( 1 1 ) 

SSD=0. 

IT0TAL=0 

IM=NCUT-1 

JM=NCUT-2 


COMPUTATION  OF  THE  COEFFICIENTS  A  1 ( I )  ,  A2(I)  AND  A3  ( I ) 

CTX 1 =  ALOG 1 0 (-ALOG ( TSTD( 1  ))) 

DO  50  1=1 , IM 
PDIF  =  PW(I)-PW(I  +  1  ) 

CTX2=AL0G10(-ALOG (TSTD (1  +  1  )) ) 

A1(I)=(PW(I ) *CTX2-PW (1+1 )*CTX1 )/PDIF 

A2(I)=(CTX1-CTX2)/PDIF 

A3 ( I  )  =  0 . 

CTX1 =CTX2 

SDTC  UT ( I )  =  0  . 

50  ITC UT ( I )  =  0 

THE  FIRST  AND  LAST  VALUES  OF  TSTD  AND  PW  ARE  CHANGED 
FOR  THE  TABLE  OUTPUT.  TRUE  VALUES  ARE  TEMPORARY  STORED 
IN  THE  RESERVE. 

TRES 1 =TSTD ( 1 ) 

TRES2=TSTD(NCUT) 

PWRES1 =PW( 1 ) 

PWRES2  =  PW ( NC  UT  ) 

TSTD ( 1 )= 1 . 0 
TSTD(NCUT )=0 . 0 
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PW ( 1 )=-1 . E70 
PW ( NC  UT  )  =  1 . E70 

PRINT  OUT  OF  THE  RESULTS 

WR ITE ( 6 , 2  )  (NAME(I) ,1=1 ,20) 

2  F0RMAT(1H1 , T25 , 20A  4 , //,T15, 'PIECEWISE -ANALYTICAL  STANDARD', 

*  'TRANSMISSION  FUNCTION  ',//,  T20 ,' TAU (X )  =', 

*  '  EXP (- 1 0** (  A 1  +  A2*X  )  )' ,///,T15, 'DATA: • ,T23, 'FROM  (  TAU 

*  ’X-VALUE)  TO  (  TAU  .X-VALUE)  WITH  (  A 1  A2  )  '  ) 

WRITE (6, 3)  (TSTD(I ) ,PW( I ) ,  TSTD ( I  + 1 ) , PW ( I  + 1 ) , A  1 ( I ) , A2 ( I ) , I  =  1 , IM ) 

3  FORMAT ( 1  HO , T28 ,  ' ( '  ,F6 . 3 , '  , '  ,F7 - 3 , ' )  (',F6.3,'  ,  '  ,  F7 . 3 ,  ' ) ' , T74 , 

*  '('  F7 . 4  '  '  F7 . 4  ')') 

WRITE (6 , 4)  ( I | WN ( I ) , CSTD ( I ) ,1=1 , NC) 

4  FORMAT(1HO,//,T15, 'ABSORPTION  BANDS: ',T40, 

*  '#  WAVENUMBER  C-VALUE ' /( 1  HO , T39 , 12 , 5X , F7 . 1 , F 1 1 . 5 ) ) 

CHECK  IF  ANY  DATA  IS  AVAILABLE  FOR  S.D.  COMPUTATION 
DATA  FORMAT  IS  THE  SAME  AS  THAT  FOR  MAIN  PROGRAMME 
ONE  CONTROL  CARD  IS  READ-IN  FIRST  FOR  BRANCHING 
IF Q  >  0  DATA  SET  FOLLOWS,  READ-IN  DATA 

IFQ  =  0  END  OF  DATA,  GO  TO  THE  FINAL  PRINTING 

READ(5 ,11) END  =  42  )  FQ , IFQ 

11  F0RMAT(5X,F10.3)T4l ,14) 

IF(IFQ.GT.O)  GO  TO  18 

42  WRITE ( 6 , 4  1  ) 

41  FORMATC///, 1H0.T5, ’$$$  NO  DATA  FOR  STANDARD  DEVIATION  COMPUTATION 

*  $$$’) 

GO  TO  40 

8  READ(5 ,11, END  =  30 )  FQ, IFQ 
IF(IFQ.LE.O)  GO  TO  30 
18  CONTINUE 
ST  =  0  . 

DO  51  1=1 , IM 
SDCUT (I ) =0 . 

51  ICUT ( I  )  =0 
CLOG= 100. 

DO  52  1=1 , NC 

52  IF( ABS(FQ-WN(I ) ) . LT . 1 . )  CLOG=CSTD(I) 

IF ( CLOG . LT . 99 . )  GO  TO  1 3 

THE  READ-IN  WAVENUMBER  DOES  NOT  MATCH  THE  MAJOR  BAND 
WAVENUMBER  (WN(I)).THE  DATA  IN  THIS  BAND  ARE  IGNORED. 

WR ITE ( 6 , 12)  FQ 

12  FORMAT( 1H0.T10, '**  ERROR  IN  WAVENUMBER  **', 

»  '  (READ-IN  WAVENUMBER  =',F10.5,'  )') 

DO  61  I DUM  =  1 , IFQ 
R E AD ( 5 , 60  )  DUMMY 

60  FORMAT ( F 1 .0) 

61  CONTINUE 
GO  TO  8 

VALID  DATA  INPUTS , READ-IN  OF  THE  DATA  AND  STAND,  RD 
DEVIATION  COMPUTATION  ARE  PERFORMED  SIMULTANEOUS, Y . 
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13  CONTINUF 
WKIliilb,  17)  FQ 

17  FORMAT ( 1 H 1 , T 1 5 , ' (  WAVE  NUMBER  =',F8.1,'  )',///, 6X , 1 WAVEN .', 3X , 

*  'PRESS. ' ,4X, 'TEMP. ' , 7X , 'U ' , 8X , 'TRANSM.  -  T(COMP)  =  DIFF’,6X, 

*  'DIFF**2 ' ,4X, 'X-VALUE  '  ,/) 

DO  9  M  = 1 , NDATA 

READ(5, 10)  EGAS, FQ, PRES, TEMP, UG,TX 
10  FORM AT (12, F10. 3, Ell . 4 ,F9 . 3 , 24X, El  1 . 4 ,F7 . 4) 

UG=UG/CF 

X  =  CLOG+AN* ALOG 10 (PRES/ 1013. )+AM*AL0G1 0(273. 1 5/TEMP ) +ALOG 1 0 ( UG ) 

DO  14  J=1 , JM 

IF(X.LE.PW(J+1 ))  GO  TO  15 

14  CONTINUE 
J  =  IM 


15  TC=EXP(-10**(A1 (J )  +  A2 ( J )  *X ) ) 


D=TX-T  C 

SD=D*D 

ST=ST+SD 

SDC UT ( J ) =SDC UT ( J  )  +SD 
ICUT ( J )  =  IC  UT ( J )  + 1 

WRITE (6, 16)  FQ,PRES,TEMP,UG,TX,TC,D,SD,X 

16  FORMAT (1H  , 3X,F8.1,F10.2,F9.2,E13.4,F9.4,F12.4,F13.6,E12.3,F9.3) 

9  CONTINUE 

END  OF  DATA  READ-IN  FOR  THIS  BAND. 

TOTAL  STANDARD  DEVIATIONS  ARE  COMPUTED  AND  PRINTED. 

20  S3D=SSD+ST 
ITOTAL=ITOTAL+IFQ 
ST=SQRT (ST/ FLOAT (IFQ) ) 

DO  21  1=1 , IM 

SDTCUT ( I ) =SDTCUT( I )+SDCUT ( I ) 

ITCUT(I)=ITCUT(I)+ICUT(I) 

21  SDCUT ( I )  =  SQRT ( SDCUT ( I ) /FLOAT ( ICUT ( I )  )  ) 

WRITE (6 ,22)  (I , TSTD ( I ) , TSTD ( I + 1 ) , ’ C UT ( I ) , SDCUT ( I ) , I = 1 , IM ) 

22  FORMAT(1H0,///,T10, 'CUTWISE  STANDARD  DEVIATION ' , // , T1 5 , ' # ' , T20 , 

*  '(  FROM,  TO  )  ',T40,'#  OF  DATA ', T53 CUTWISE  SD ' , // , ( T 1 4 , 12 , 

*  T20, ' (  '  ,F5.2, ' , ' ,F5.2, ' ) ' ,T43, I  4 , T52 , F 1 0 . 6 , / )  ) 

WRITE ( 6 , 23 )  IFQ, ST 

23  FORMAT(  1 H  0  ,  T 1  0  ,  'TOTAL  it  OF  DATA  FOR  THIS  BAND  =',I5,9X, 

*  '  STANDARD  DEVIATION  =',F12.6,//) 

GO  TO  8 

END  OF  THE  STANDARD  DEVIATION  COMPUTATION  FOR  ALL  DATA. 

GRAND  TOTAL  STANDARD  DEVIATION  IS  COMPUTED  AND  PRINTED  OUT 
TOGETHER  WITH  VITAL  INFORMATIONS. 

30  SSD=SQRT(SSD/FLOAT(ITOTAL) ) 

DO  31  1  =  1  ,  IM 

31  SDTCUT (I )=SQRT (SDTCUT (I ) /FLOAT ( ITC UT ( I ) )  ) 

WRITE (6, 32)  (I , TSTD(I ) , PW( I ) ,TSTD(I  +  1 ) , PW ( I  + 1 ) . A  1 ( I  )  , A2 ( I  )  , 

*  ITCUT(I) ,SDTCUT(I)  ,1  =  1 ,IM) 

32  FORMAT( 1H  1  ,T20, '**»  P IECEWISE-ANALYTICAL  STANDARD  TRANSMISSION', 
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*  '  FUNCTION  ***',///, T10, 'TOTAL  CUTWISE  STANDARD  DEVIATION',//, 

*  T  15  ,  'CURVE  // '  ,  3X  ,  '  FROM  (  TAL  .X-VALUE)  TO  (  TAU  .X-VALUE)', 

*  '  WITH  (  A 1  A2  )  ' , 3X, ' #  OF  DATA ', 4X ,' CUTWISE  SD',//, 

*  (T18,I2,T30,  '  (  '  ,F6.3,  '  ,’,F7.3,'>  C,F6.3,’  ,  '  ,  F7 . 3  ,  '  )  '  ,  T76  , 

*  ’  (  '  ,F7.4, ' , ' ,F7.4, ' ) ' ,7X, 1 3, 5X, FI 0.6,/) ) 

WRITE (6 , 33 )  ITOTAL , SSD 

33  FORMAT(1H0,T10, 'GLOBAL  RESULTS ',//, T 1 5 TOTAL  NUMBER  OF  DATA', 

*  '  USED' ,15,//, T15, 'GLOBAL  STANDARD  DEVIATION' ,F12. 6) 

40  CONTINUE 

END  OF  ALL  COMPUTATION. 

RESERVED  TRUE  VALUES  OF  THE  FIRST  fl ND  LAST  TSTD 
AND  PW  ARE  RETURNED. 

TSTD ( 1 ) =TRES 1 
TSTD(NCUT)=TRES2 
PW( 1 ) =PWRES 1 
PW(NCUT)=PWRES2 

CALL  SDTAU 

RETURN 

END 

SUBROUTINE  INTPL2 

COMPUTATION  OF  THE  PIECEWISE-ANALYTICAL  STANDARD  TRANSMISSION 
FUNCTION 

VERSION  2  -  1 

TAU  =  EXP(-10**(A1+A2*X+A3*X**2‘,  ) 

COMMON  /PARM 1 /  TSTD ( 1 2 ) , PW ( 1 2 ) , WN ( 6 ) , CSTD ( 6 ) , NC UT , NC , NAME ( 2 0 )  , 

»  AN, AM,CF, ICONST(6) ,NEL 

COMMON  /PARM3/  A  1 ( 1 1 ) , A2 ( 1 1  )  , A3 ( 1 1  ) 

DIMENSION  T( 10) ,TD(6 ,74) ,PD(6 ,74)  ,  JI (6 , 10)  , 

*  SDK (7) , SDE ( 7 , 9 ) ,SUME(2,9) ,DE(6,9) 

NWC  =  0 

K  =  0 

READ(5 , 2 , END=80 )  FREQ , MAXDAT 
IF ( MAXDAT . GT . 0 )  GO  TO  3 
80  CONTINUE 
WRITE (6 , 99  ) 

99  FORMAT (/// , 1  HO , T5 ,  '$$$  NO  DATA  FOR  STANDARD  DEVIATION  COMPUTATION 

*  $$$') 

GO  TO  77 

1  CONTINUE 

READ(5,2,END  =  21  )  FREQ, MAXDAT 

2  FORMAT (5X , F 1 0 . 5 , T4 1 ,14) 

IF ( MAXDAT . EQ . 0  )  GO  TO  21 

3  CLOG= 1 . E  10 
DO  5  L=1 , NC 

IF ( ABS ( FREQ-WN (L ) ) . LE .  .01)  CLOG=CSTD(L) 

5  CONTINUE 

IF ( ABS ( CLOG- 1.E10).GE..01)  GO  TO  9 
WR ITE (6 , 100)  FREQ 
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100  FORMAT ( '  1 ’ ,///////, '  ERROR  IN  INPUT  DATA  ;  WAVE  NUMBER  ’  , F 1 0 . 3 , 

*  '  NOT  USED  IN  COMPUTATION  OF  CONSTANTS.') 

DO  6  J  =  1  ,  MAXDAT 

READ(5 , 101 )  KGAS, FREQ, PRES, TEMP, PPM, RANGE, UGAS.TX 
6  CONTINUE 
GO  TO  1 
9  CONTINUE 
K=K  + 1 
NWC=NWC  + 1 
J I ( K , 1  )  =  0 
JI (K , NCUT )=MAXDAT 
J  =  2 

DO  20  1=1 , MAXDAT 

READ(5, 101)  KGAS, FREQ, PRES, TEMP, PPM, RANGE, UGAS.TX 

101  FORMAT(I2,F10.3,E11 .4,F9.3,E11.4,E13.6,E11.4,F7.4) 

TD (K , I ) =TX 

PD(K , I )=  AN  *AL0G10(PRES/1013- )+  AM  *ALOG 1 0 (273 . 1 5/TEMP )+ALOG 1 0 

*  ( UGAS/CF ) +CLOG 
IF(TD(K,I).GE. TSTD( J ) )  GO  TO  20 
IF(J.EQ.NCUT)  GO  TO  20 
JI(K,J)=I— 1 

J  =  J  +  1 

20  CONTINUE 
GO  TO  1 

21  CONTINUE 
IF(NWC.LE.O)  RETURN 
DO  30  J  =  1 , NC UT 

T(J )=ALOG10(-ALOG(TSTD(J )  )  ) 

30  CONTINUE 
SUMT  =  0  . 

DT  =  0 . 0 
NC  C  =  NC  UT - 1 
DO  45  1=1 , NCC 
SA=0  . 

TA=0  . 

UA=0 . 

DO  41  K= 1 , NWC 
SUME(K,I)=0.0 
M=JI (K , I )+1 
N  =  J I (K , 1  + 1  ) 

DO  40  J  =M  ,  N 

TC  =  ALOG 10 ( -ALOG (TD (K , J  ) )  ) 

SA=SA+  TC  *PW(I)*PW(I+1 )*(PD(K, J )-PW( I ) ) * (PD (K , J)-PW(I  +  1 )  ) 

TA=TA+( (PD(K, J)-PW(I) )*(PD(K, J)-PW(I+1 ) ) * ( (PD (K , J ) **2 ) « (PW( I + 1 ) 

*  *T(I )-PW(I )*T(I+1 ) )+PD(K , J)*( (PW(I )«*2)*T(I+1 )-(PW(I+1 )**2)»T(I) 

*  ) )  )  / ( PW( I )-PW( 1  +  1 ) ) 

UA=UA+((PD(K,J)-PW(I))*(PD(K,J)-PW(I+1)))**2 

40  CONTINUE 

41  CONTINUE 

A  1 ( I )  =  ( SA-T  A  )  /U A 

A3(I)=(T(I)-T(I+1 ))/(PW(I+1 )*(PW(I)-PW(I+1 )))-(T(I)-A1(I))/(PW(I)* 

*  PW(I+1 ) ) 

A2(I  )= (T (I )-T (1  +  1 ) )/ (PW( I ) -PW ( I  + 1 ) )-A3(I )* ( PW( I )+PW(I  +  1 ) ) 

D I  =0  . 

SUM  1  =  0  . 

DO  44  K  =  1  ,  NWC 


N=JI (K, 1+1 ) 

DE (K , I ) = FLOAT ( 1 +N-M) 

DO  43  J=M,N 

SUME(K, I)=SUME(K,I)+(TD(K, J)-EXP(-10.**(A3(I)*PD(K, J)*PD(K, J)+ 

*  A2(I)*PD(K,J)+A1(I))))**2 

43  CONTINUE 

SDE(K,I)=SQRT(SUME(K,I)/DE(K,I)) 

SUM I=SUMI+SUME (K , I )* FLOAT ( ICONST (K  ) ) 

DI =DI+DE (K , I ) * FLO AT ( ICONST (K ) ) 

44  CONTINUE 
SUMT=SUMT+SUMI 
DT=DT+DI 

SDE (NWC  + 1 , I)=SQRT(SUMI/DI ) 

45  CONTINUE 

DO  51  K= 1 , NWC 
SUMK=0 . 0 
DK=0 . 0 

DO  50  1=1 , NCC 
SUMK=SUMK+SUME(K, I) 

DK=DK+DE (K , I ) 

50  CONTINUE 
SDK(K)=SQRT(SUMK/DK) 

51  CONTINUE 

SDKCNWC+1 ) =SQRT (SUMT/DT ) 

DUM 1 =PW ( 1  ) 

DUM2  =PW ( NC  UT ) 

DUM3=TSTD( 1  ) 

DUM4=TSTD(NCUT) 

PW( 1 )=-1000000. 

PW(NCUT)= 1000000. 

TSTD ( 1 )=1 .0 
TSTD(NCUT )  =  0 . 0 
DO  60  K= 1 , NWC 

WRITE (6, 1 02 ) ( NAME ( J ) , J  =  1 , 20 )  , WN (K ) , NCC , ( I , TSTD ( I ) , TSTD ( 1  +  1 ) ,PW(I) , 

*  PW(I+1 ) ,A1(I) ,A2(I) , A3 ( I ) ,SDE(K,I) ,1=1 ,  NCC) 

102  FORMAT  (' 1 35X RENDITION  OF  EMPIRICAL  TRANSMITTANCE  FUNCTION 
*FOR  :  '//,20A4,//,40X, 'WAVE  NUMBER  : '  , F 1 5 . 4 , //// , 20X , ' THE  TRANSMI 
*SSION  CURVE  IS  DIVIDED  INTO', 13,'  SEPARATE  CURVES /20X EACH  CU 
*RVE  IS  EXPRESSED  BY  A  FUNCTION  OF  THE  FORM  "  TAU  =  EXP ( - 1 0*« ( A3*P 
//*P+A2*P+A  1  ) )  "  .  '  /  ,  20X  ,  'THE  FUNCTION  COEFFI 

•CIENTS  AND  RESULTING  STANDARD  DEVIATION  FOR  EACH  CURVE  ARE  AS  FOL 
*LOWS : ' ,////22X, ' TAU ' , 20X , ’ P ' ,24X, '  A  1  '  ,  1  3X, ' A2 ' , 1 3X, ' A3 ' , 1 1X, 'STAND 

*  ARD  DEVIATION',///,  (5X, 'CURVE  # ' , 13, 3X , • ( ' ,F4. 2, ' - '  ,F4. 2, ' )  '  , 5X, 
•'(' ,F9.5, '-' ,F9.5, ') ' ,5X,3F15.6,  6X,F15.6  /)) 

WR ITE ( 6 , 40 1  )  SDK (K ) 

401  FORMATdX,// ,87X, 'TOTAL  STANDARD  DE VI ATION ’ , F 1 5 . 6 ) 

60  CONTINUE 
K  =  K  + 1 

WRITE ( 6 , 104)  (NAMECJ ) , J=1 ,20) ,NCC,  ( I , TSTD ( I ) , TSTD ( 1+ 1 ) , PW( I ) , 

*  PW ( 1  + 1 ) , A  1 (I ) , A 2 ( I ) , A 3 ( I ) , SDE (K , I) , 1  =  1 , NCC) 

104  FORMAT  C 1 ',/, 35X ,' RENDITION  OF  EMPIRICAL  TRANSMITTANCE  FUNCTION 
*FOR  :  ’//,20A4,//,40X, 'TOTAL  PROFILE  AVERAGED  OVER  ALL  WAVE  NUMBER 
*S'  ,////, 20X, 'THE  TRANSMI 

*SSION  CURVE  IS  DIVIDED  INTO’, 13,'  SEPARATE  CURVES .', /20X ,' EACH  CU 
•RVE  IS  EXPRESSED  BY  A  FUNCTION  OF  THE  FORM  "  TAU  =  EXP ( - 1 0*« ( A3*P 
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#*P+A2*P+A 1 ) )  " . ' / , 20X ,  'THE  FUNCTION  COEFFI 

•CIENTS  AND  RESULTING  STANDARD  DEVIATION  FOR  EACH  CURVE  ARE  AS  FOL 
•LOWS:  ' ,////22X, 'TAU ' , 20X , ' P • , 24X , ' A  1 ' , 1 3X, ’ A2 ' , 1 3X , ' A3  '  , 1 IX , 'STAND 
*ARD  DEVIATION',///,  (5X, 'CURVE  # * , 13 , 3X , • ( ' ,F4 . 2 , ' - ' ,F4 . 2 , ' ) ' , 5X , 
*'( '  ,F9.5, ,F9.5, ' ) ' , 5X , 3F 1 5 . 6 ,  6X.F15.6  /)) 

WRITE (6 , 402  )  SDK (K ) 

402  FORMAT(1X,//,81X, 'GRAND  TOTAL  STANDARD  DEVIATION ', F 1 5 . 6 ) 

WRITE (7,201 ) ( A 1 (I) ,A2(I) , A  3 ( I ) ,1  =  1 ,NCC) 

201  FORMAT ( 3F 10.6) 

IDT=IFIX(DT ) 

WRITE (6, 225) IDT, NEL, SDK (K) 

225  F0RMAT(// ,  1H0,T15  ,  'GRAND  TOTAL  //  OF  DATA  =  '  ,  15  ,  /  /  ,  T 1 5  ,  '  #  OF', 

*  '  ELIMINATED  DATA  = ' , 1 5 , // , T 1 5 , 'GLOBAL  STANDARD  DEVIATION  IN', 

*  '  TAU  =' , F 12 . 6 , // ) 

DO  76  M  =  1 , NWC 

IF(ICONST(M) .EQ.1)  GO  TO  76 
WRITE (6 , 226  )  M 

226  FORMAT ( 1 H  ,T15,'N0TE:  THE  BAND', 13,’  IS  NOT  INCLUDED  IN  THE', 

*  '  FINAL  STANDARD  DEVIATION') 

76  CONTINUE 
PW(1 ) =DUM 1 
PW ( NC  UT )  =  DUM2 
TSTD ( 1 )=DUM3 
TSTD (NC  UT ) =DUM4 

CALL  SDTAU 

77  CONTINUE 
RETURN 
END 

SUBROUTINE  SDTAU 

COMPUTATIONS  OF  STANDARD  DEVIATIONS  IN  TAU  USING  THE  ORIGINAL 
DATA  USED  IN  MAIN 

DIMENSION  NDATA( 6 ) ,TSD(6) ,WWW(12, 10) ,STANDV(12) ,P(10),T(10) 

COMMON  /PARM 1 /  TSTD ( 1 2 ) , PW ( 1 2 ) , WN ( 6 ) , CSTD (6 ) , NC UT , NC , NAME (20 ) , 

*  AN,AM,CF,ICONST(6) ,NEL 
COMMON  /PARM2/  PRES (6 , 1 2 , 1 0 )  , TEMP (6 , 1 2 , 1 0 ) , UGAS (6 , 1 2 , 1 0 ) , 

*  TAU(6, 12) , NTC ( 6 ) , NLV (6 ) 

COMMON  /PARM3/  A 1 ( 1 1 ) , A2 ( 1 1 ) , A3 ( 1 1 ) 

C 

NGDATA  =  0 
GTSD=0. 

ICST=NC 
DO  70  M= 1 , NC 
JM  =NTC ( M) 

KM=NLV( M) 

NDATA ( M) = JM*KM 

NGDATA =NGD AT A +N DATA (M)*IC0NST(M) 

TSD( M) =0 . 

WRITE(6 ,214) 

214  FORMATC1  ',/////, 45X, 'RECOMPUTATION  OF  TAU',/////) 

IF(M.GT. 1 )  GO  TO  77 
WRITE (6 ,215) 


o  n  o  o  o  o 


215  FORMAT (20X , 'A  TAU  VALUE  ,  T  ,  IS  RECOMPUTED  FOR  THE  ORIGIONAL  DATA 

*  USING  THE  PIECEWISE -A  NAUTICAL  TRANSMISSION  FUNCTION  //20X , 

*  'STANDARD  DEVIATIONS  BETWEEN  THE  ACTUAL  TAU  AND  THE  RECOMPUTED’, 

*  ’  TAU  VALUES  ARE  COMPUTED.’////) 

77  CONTINUE 

WRITE (6 , 202 )  M,WN(M) ,NTC(M) ,NLV(M) 

202  FORMAT(1HO,T15, ’**#  CASE’, 13,’  (WAVE  NUMBER  =’,F10.3,’)  •••’, 

*  ///,T20, ’TOTAL  #  OF  CUTS  =’, 13 ,//, T20 ,’ TOTAL  #  OF  LEVELS  =’,I3, 

*  ///) 

WRITE (6 ,216)  AN , AM , M , CSTD( M) 

2 1 6  FORMAT ( 1  OX , ’ N  =  ’  , F 1 0 . 5 , // 1  OX , ’ M  = ’ , F 1 0 . 5 , / / , 1 0X , ’ C ’ , 1 1 , ’  =’, 

*  FI 0 . 5 , ////) 

WRITE (6 ,217)  - 

217  FORMAT(//,1HO,T7, ’RECOMPUTED  TAU  AND  STANDARD  DEVIATIONS’, 

*  ’  IN  TAU  ’ ,/, 1H0,T2, ’CUT* ,T11 , ’TAU’ ,T20, ’X*’ ,T30, ’XI ’ ,T39, 

*  ’X2’ ,T48, ’X3 ’ ,T57, ’X4 ’ ,T66, ’X5 ’ ,T75, ’X6 ’ ,T84, ’X7’ ,T93, ’X8’ , 

*  T102, ’ X 9 ’ ,T111 , *  X 1 0 ' , T 1 2 1 , ’ CUTWISE-SD • ,/) 

COMPUTATION  OF  THE  CUTWISE  STANDARD  DEVIATIONS  IN  X 

DO  71  J=1 , JM 
WW=0. 

DO  72  K= 1 , KM 

P(K)=CSTD(M)+AN#PRES(M,J,K)+AM*TEMP(M,J,K)+UGAS(M,J,K) 

IM=JM-1 
DO  75  1=1 , IM 

IF(PW(I+1 ) .GT.P(K) )  GO  TO  76 

75  CONTINUE 
I  =  IM 

76  CONTINUE 

T(K)=EXP(-10»«(A3(I)*P(K)«P(K)+A2(I)*P(K)+A1 (I))) 

WWW( J,K)=(TAU(M, J)-T(K))»*2 
WW=WW+WWW( J ,K) 

72  CONTINUE 

WW=SQRT ( WW/FLOAT (KM ) ) 

WRITE (6, 21 8)  J , PW ( J ) , T A U ( M , J ) ,(T(K)  ,K=1 ,KM) 

218  FORMAT ( 1H  , 15 , F9 . 4 ,F9 . 4, IX , 10F9 . 4) 

WRITE (6 ,219)  WW 

219  FORMAT(1H+,T121 , F 1 0 . 5 ) 

71  CONTINUE 

COMPUTATION  OF  THE  LEVELWISE  STANDARD  DEVIATIONS  IN  X 

DO  73  K= 1 , KM 
WW=0. 

DO  74  J=1 , JM 
WW=WW+WWW ( J , K ) 

74  CONTINUE 

TSD(M)=TSD( M) +WW 

STANDV(K)  =  SQRT (WW/FLOAT ( JM)  ) 

73  CONTINUE 

WRITE (6 , 220 )  (STANDV(K) ,K=1 ,KM) 

220  FORMAT( 1  HO , T4 , ’ LE VELWISE-SD  : ’ ,T26 , 10F9. 5) 

GTSD=GTSD+TSD( M) * FLOAT ( ICONST ( M) ) 

ICST=ICST-ICONST(M) 

TSD(M)=SQRT(TSD(M) /FLOAT (NDATA(M) ) ) 
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WRITE (6,221)  TSD(M) 

221  F0RMAT(1H0,//,T15, 'TOTAL  STANDARD  DEVIATION  FOR  THIS  CASE 

*  F15.6) 

70  CONTINUE 

GTSD=SQRT (GTS D/ FLOAT (NGDATA  ) ) 

WRITE (6 , 223 )  AN, AM 

223  FORMAT ( ’  1 ' , T15 , ’  ***  SUMMARY  OF  THE  TRANSMITTANCE  RECOMPUTATION  * 

*  **’,/// ,T20, ’PRESSURE  EXPONENT  N  = • , F 1 0 . 5 , // , T20 , 

*  ’TEMPERATURE  EXPONENT  M  = ’  ,  F 1 0  5 , /// , T5 , ’ CASE  #’,3X, 

*  ’WAVE  NUMBER’ ,5X, ’C-VALUE  ’  ,5X, ’TOTAL  //  OF  DATA’,3X, 

»  ’CA^FWT^F  n  TN  TAIMI 

WRITE (6 ,224) ' (M,WN(M) ,CSTD(M) , NDATA( M) , TSD(M) , M= 1 , NC ) 

224  FORMAT(1HO,T6,I3,6X,F9.2,5X,F8.3, 10X.I3, 12X.F12.6) 

WRITE (6 , 225  )  NGDATA , NEL , GTSD 

225  FORMATC// ,  1H0,T15 , ’GRAND  TOTAL  #  OF  DATA  = ’ , 15, // ,T15 , ’ #  OF’, 

*  ’  ELIMINATED  DATA  =’, 15 ,//, T 15 ,’ GLOBAL  STANDARD  DEVIATION  IN’, 

*  ’  TAU  =’  ,F12 . 6  ,  //) 

IF(ICST.LE.O)  RETURN 
DO  78  M=1 , NC 

IF ( ICONST(M) . EQ. 1 )  GO  TO  78 
WRITE (6 , 226 )  M 

226  FORMAT ( 1 H  ,T15,’NOTE:  THE  BAND’, 13,’  IS  NOT  INCLUDED  IN  THE’, 

*  ’  FINAL  STANDARD  DEVIATION’) 

78  CONTINUE 

RETURN 

END 


174 


C 

C**» 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


***  COMPUTER  CODE  SIMMIN  *»»*** 

VERSION  (  6  -  3  )  TRACE  GASSES 

COMPUTATION  OF  ABSORBER  PARAMETERS  AND  ANALYTICAL  STANDARD 
TRANSMISSION  FUNCTION 

THIS  CODE  USES  THE  SUBROUTINE  FMCG  IN  SSP  LIBRARY 
THIS  CODE  CONSISTS  OF 

MAIN:  DATA  READ-IN  AND  CONTROL  OF  COMPUTATION 

FUNCT:  COMPUTATION  OF  THE  COST  AND  ITS  DERIVATIVES 

TITLE:  PRINTING  OF  HEADINGS  AND  INITIAL  CONDITIONS 

PRINT1 :  PRINTOUT  OF  RESULTS  AND  COMPUTATION/PRINTING  OF  S.D.S 

NMBC :  COMPUTATION  OF  NON-MAJOR  BANDS'  C-VALUES 


DATA 
1  . 


2. 


3- 

4. 


SET-UP 

INITIAL  GUESSES  X(I)  (9  CARDS  WITH  T12.F10.7) 

X ( 1 ) =A 1 ,  X ( 2 ) = A2 ,  X(3)=A3,  X(4)=N,  X(5)=M,  X(5+I )=LOG(C( I ) ) 
(NEED  DUMMY  INPUTS  FOR  PROBLEMS  WITH  DIMENSION  <  9) 

SIGNAL  VARIABLES  S ( I )  (9  CARDS  WITH  T12.F10.7) 

SCI)  =  0  ...  X ( I )  IS  KEPT  CONSTANT 
SCI)  =  1  ...  X ( I )  IS  VARIED 

(NEED  DUMMY  INPUTS  FOR  PROBLEMS  WITH  DIMENSION  <  9) 

COMMENT  CARD  (20A4)  FOR  TITLE  AND  ABSORBER  TYPE  ETC. 

ONE  FOR  EACH  ABSORPTION  BAND 


#  OF  DATA  AND  COMMENTS 
101  ) 

102) 


DATA  SETS  (MAX.  4  SETS)  - 
EACH  SET  CONSISTS  OF 

1  ST ( CONTROL )  CARD:  WAVENUMBER, 

(SEE  FORMAT 

DATA  CARDS:  P,  T,  U,  TAU  ETC. 

(SEE  FORMAT 
(TOTAL  #  OF  DATA  SHOULD  NOT  EXCEED  900) 

5.  BLANK  CARD  -  FOR  THE  TERMINATION  OF  DATA  INPUT  FOR  MAIN 

6.  DATA  SETS  FOR  NMBC  -  ONE  FOR  EACH  ABSORPTION  BAND 

EACH  SET  CONSISTS  OF 

DATA  CARDS:  SAME  AS  MAIN 
FINAL  CARD:  BLANK 

TERMINATION  A  CONTROL  CARD  WITH  -1  IN  FIRST  TWO  COLUMNS 
THIS  COMES  AFTER  THE  FINAL  BLANK  CARD 
(IF  NO  DATA  BUT  A  BLANK  CARD  IS  SUPPLIED,  NMBC  IS  SKIPPED) 

DIMENSION  X(9),G(9),Y(9),H(72),WN(4) 

COMMON  /PARM 1 /  NC , ND(5 ) , RW( 4 ) 

COMMON  /PARM2/  IC , PLOG ( 900 ) , TLOG ( 900 ) , ULOG ( 900 ) , TAU ( 900 ) , S( 9 ) 
COMMON  /PARM3/  P ( 900  )  ,  T ( 900 )  , U ( 900 ) , L (20 ) 

COMMON  /PARM4/  P 0 , TO , NDIM , ID ( 5 , 9 ) 

EXTERNAL  FUNCT 

CONSTANTS 
P0=1 . 0 1 3E+03 
T0=273 • 15 
CF=2 . 69E+ 1 9 
N  =  9 
V  =  0. 

IC  =  0 
MAXNC=4 


bJ 
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DATA  INPUT 

READ (5 , 100)  (X(I),I=1,N) 

READ(5 ,100)  (S(I) ,1=1 ,N) 

100  FORMAT  (T12.F10.7) 

COMMENT  CARD  (THIS  INCLUDES  THE  ABSORBER  TYPE) 

READ(5 , 500 )  (L ( I ) , 1= 1 , 20 ) 

500  FORMAT (20A4 ) 

NC  =  #  OF  MAJOR  ABSORPTION  BANDS 

ND ( 1 ) =0 ,  ND ( 2 ) =N  1  ,  ND ( 3 )  =  N 1 +N2 ,  ND(4 )  =  N 1+N2+N3,  ... 

WHERE  N 1  ,  N2 ,  N3,  ...  ARE  //S  OF  DATA  IN  BANDS  1,  2,  3,  . 
ND(NC  +  1  )=  TOTAL  0  OF  DATA 

NC  =  0 
ND ( 1  )  =0 

DO  10  M= 1 , MAXNC 

READ (5 , 101 )  WN(M) ,IX, (ID(NC+1 , I) , 1=1 ,9) 

101  FORMAT(5X,F10.3,T41 ,I4,9A4) 

IF(IX.LE.O)  GO  TO  11 
NC=NC+1 

IM=ND (NC ) + 1 
IN=ND(NC )+IX 
ND(NC+1 )=IN 
DO  12  I  =IM , IN 

READ (5, 102)  P (I ) ,T(I) ,U(I) ,TAU(I) 

102  FORMAT  ( 1 2X , E 1 1 . 4 , F9 . 3 , 24X , E 1 1 . 4 , F7 . 4 ) 

U(I)=U(I)/CF 

DATA  ARE  CONVERTED  TO  THE  LOG  OF  THE  NORMALIZED  VALUES 

PLOG(I)=ALOG10(P(I)/PO) 

TLOG ( I ) = ALOG 10(T0/T(I)) 

ULOG ( I ) = ALOG 1 0 (U ( I ) ) 

12  CONTINUE 

10  CONTINUE 

11  CONTINUE 

END  OF  DATA  INPUT 

IF(NC.GT.O)  GO  TO  13 
WRI TE ( 6 , 110) 

110  FORMAT  (1H0, 'ERROR  IN  DATA  INPUT') 

GO  TO  1000 

13  CONTINUE 

NDIM=0 
N  =5+NC 
DO  14  1=1 , N 

IF ( S ( I )  .  NE . 0 . )  NDIM  =  NDIM+ 1 

14  CONTINUE 

DO  15  1  =  1  , NC 

RW(I)=FLOAT(ND(NC  +  1 ) ) / ( FLOAT (N D ( I  + 1 ) — N  D ( I ) )*FLOAT(NC) ) 
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RW ( I ) = 1 .0 

15  CONTINUE 
C 

DO  16  1  =  1  ,N 

Y  ( I )  =  X ( I ) 

16  CONTINUE 
C 

EST  = 1 . E-6 
EPS= 1 . E-6 
LIMIT=1 
I  ER  =  0 
C 

c#«»***»*«*##*  fmcg  SEARCH  a**********************  START  »**»***«**•***»» 

C 

CALL  FMCG ( FUNCT, N,X,V,G, EST, EPS, LIMIT ,IER,H) 

C 

c »****»»*«*»»»  fmcg  SEARCH  ***********************  END  **************** 

C 

CALL  TITLE(N,Y, LIMIT, EPS) 

C 

CALL  PRINT1 (N ,X, V,G, IER) 

C 

CALL  NMBC(X,L,NC,WN,CF,PO,TO) 

C 

1000  CONTINUE 
STOP 
END 

SUBROUTINE  FUNCT ( N , X , V , G ) 

COMPUTATION  OF  THE  FUNCTION  VALUE  AND  DERIVATIVES 

(DOUBLE  EXPONENTIAL  FUNCTION) 

DIMENSION  X ( 9 )  ,G(9)  ,F(9) 

COMMON  /P ARM  1 /  NC , ND (5 ) , RW( 4 ) 

COMMON  /PARM2/  I C , PLOG ( 900 ) , TLOG ( 900 ) , ULOG ( 900 ) , TAU ( 900 ) , S ( 9  ) 

C 

IC=IC+1 

V  =  0. 

DO  20  K= 1  ,  N 
G (K ) =0 . 

20  CONTINUE 
C 

DO  21  J=1 ,NC 
J J= J+5 
SQE  R  =  0 . 

DO  22  L  =  1  ,5 
F ( L ) =0 . 

22  CONTINUE 
F  (  J J  )  =0 . 

IM=ND(J )+1 
IN=ND( J+1 ) 

DO  23  I =IM , IN 

W1=X( JJ  )+X(4)*PL0G(I )+X(5 )*TLOG(I )+ULOG(I ) 

R=X( 1 )+X(2)«W1+X(3)*W1*W1 
R  =  1 0 . *  *R 
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IFCR.LE.70. )  GO  TO  24 
TC  =  0  . 

GO  TO  25 

24  CONTINUE 
TC  =  EXP(-R  ) 

25  CONTINUE 
E=TAU(I)-TC 
R=R*E*TC 

SQE  R  =  SQER+E  **2 
F ( 1 )=F( 1 )+R 
F(2)=F(2)+R*W1 
F( 3 )=F( 3 )+R*W1 *W1 
R=R*(X(2)+2.*X(3)*W1 ) 

F(4)=F(4)+R*PLOG(I) 

F(5)=F(5)+R*TLOG(I) 

F(JJ)=F( JJ)+R 
23  CONTINUE 

V=V+SQER*RW ( J  ) 

DO  26  K=1 ,5 
G(K)=G(K)+F(K)*RW(J ) 

26  CONTINUE 
G(JJ)=F(JJ)*RW(J) 

21  CONTINUE 

DO  27  1=1 ,N 

G(I)=4. 60517 *G(I)*S(I) 

27  CONTINUE 

RETURN 

END 

SUBROUTINE  TITLE ( N , X , LIMIT , EPS ) 

PRINTING  OF  THE  TITLE  AND  INITIAL  VALUES 

DIMENSION  X ( 9 ) , L ( 4 ) 

COMMON  /P ARM  1 /  NC , ND (5  )  ,  RW ( 4  ) 

COMMON  /PARM4/  PO , TO , NDIM  ,  ID (5 , 9 ) 

DO  40  1=1 , NC 
L(I)=ND(I  +  1 ) -N  D ( I ) 

40  CONTINUE 

CALL  DATE  ( MONTH , IDAY , IYEAR  ) 

WRITE (6, 1 1 1 ) MONTH, IDAY, IYEAR 
111  F0RMAT(1H1 ,T60,I4,  '  /  ',12,'  /  ’,12,/) 

WRITE (6 , 400  )  NC , NDIM 

400  FORMAT  (1H  ,T14,'**»  SIMULTANEOUS  PARAMETER  EVALUATION  ***',///, 

*  '  PARAMETERS  :  (  N  ,  M  ,  A 1  ,  A2  ,  A3  ,  C(I),  1  =  1, ',12,'  )’, 

*  8X,'(  DIMENSION  = ' ,13, '  )',///,'  DATA  :') 

WRITE (6, 401)  ((ID(K,J),J=1,9),L(K),K=1,NC1 

401  FORMAT(  1H  +  ,T1  1  ,  '  (  '  ,9A4,  '  )’,5X,'//  OF  POINTS  =  '  ,15,//) 

WRITE (6, 402)  N D ( NC + 1 ) , PO , TO , LIM IT , EPS 

402  FORMAT  (  1  H+  ,  T5  1  ,’  TOTAL  if  OF  DATA  =',I5,///, 

*  '  FUNCTION  :  TAU  (  W  )  =  EXP  (  -10  **  (  A 1  +  A2  *  W  + ' , 

*  '  A3  *  W* *2  +  A4  *  W**3  )  )',// , T 1 5 , 'WHERE,  ', 
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*  '  W  =  LOG ( C )  +  LOG (  U  *  (P/P0)**N  «  (T0/T)»*M  ) 

*  '  ,  A4  =  0  .*,///  , 

*  '  CONSTANTS  :  PO  =  '  ,  F8 . 2 , 7X , ' TO  =',F8.2,7X, 

*  'LIMIT  =•  ,I5,7X, 'EPS  =  '  ,  1  PE  1 0 . 1 , / ) 

IF(NC.EQ.I)  WRITE(6,403)  X ( 1 ) , X ( 6 ) , X (2 ) , X ( 3 ) , X ( 4 ) , X (5 ) 

IF ( NC . EQ . 2 )  WRITE (6 , 404 )  X ( 1 ) , X (6 ) , RW( 1 ) , X (2 ) , X(7 ) , RW( 2 ) , X ( 3 ) , 

*  X  ( 4  )  ,  X  ( 5  ) 

IF (NC . EQ . 3 )  WRITE (6 , 405 )  X ( 1  ) , X (6 ) , RW( 1 ) , X (2 ) , X (7 ) , RW( 2 ) , X ( 3 )  , 

*  X(8),RW(3),X(4),X(5) 

IF ( NC . EQ . 4 )  WRITE ( 6 , 406 )  X ( 1 ) , X (6 ) , RW( 1 ) , X (2 ) , X ( 7 ) , RW( 2 ) , X ( 3 ) , 

*  X ( 8 )  ,RW(3)  ,X(4) ,X(9) ,RW(4) ,  X  ( 5  ) 

403  FORMAT  ( 1  HO INITIAL  VALUES  :',T22,'A1  = *  ,  F 1 2 . 7 , 9X , • LOG ( C 1 )  =’, 

«  F 1 2 . 7 , //  , T22 , ' A2  = ' , F 1 2 . 7 , // , T22 , ’ A  3  = '  ,  F 1 2 . 7 , / / , T22 , ' N  =', 

*  FI  2 . 7 , // , T22 , 'M  ='  , F 1 2 . 7  , / ) 

404  FORMAT  ( 1  HO  INITIAL  VALUES  :’,T22,’A1  =  ’  ,  F 1 2 . 7 , 9X  ,  ’  LOG  ( C 1  )  =\ 

*  F 1 2 . 7 , 4X , ' (  WEIGHT  =’,F8.4,'  )',//,  T22 ,’ A2  = • ,F 12 . 7 , 9X , 

*  'LOG ( C2 )  =’  ,F12.7,4X, ' (  WEIGHT  =’,F8.4,'  )',//, T22 A3  =', 

*  F12.7.//.T22, 'N  = ' , F 1 2 . 7 , // , T22 , ’ M  =',F12.7,/) 

405  FORMAT  ( 1H0 ,' INITIAL  VALUES  :,,T22,'A1  = ’ , F 1 2 . 7 , 9X , ' LOG ( C 1 )  =’, 

*  F 1 2 . 7 , 4X , ' (  WEIGHT  =’,F8.4,'  )',//, T22 ,' A2  =',F12.7,9X, 

*  ' LOG ( C2 )  =' ,F12.7,4X, ’(  WEIGHT  =',F8.4,’  )’,//, T22 ,' A3  =  ’, 

*  FI 2 . 7 , 9X , ' LOG ( C 3 )  = ' ,F 1 2 . 7 , 4X , ' (  WEIGHT  =’,F8.4,'  )',//, T22, 

*  'N  ='  F12.7  //  T22  'M  ='  F12.7  /) 

406  FORMAT  ( 1  HO INITIAL ’ VALUES ’  :',T22,'A1  = ' , F 1 2 . 7 , 9X , ' LOG ( C 1 )  =  ', 

*  F 1 2 . 7 , 4X , ' (  WEIGHT  =’,F8.4,'  )',//, T22 ,' A2  =',F12.7,9X, 

*  ' LOG ( C2 )  =' ,F12.7,4X, '(  WEIGHT  =',F8.4,'  )',//, T22 ,' A3  =’, 

*  F 1 2 . 7 , 9X , ' LOG ( C3 )  =  '  , F 1 2 . 7 , 4X  ,  '  (  WEIGHT  =',F8.4,'  )',//, T22, 

*  'N  ='  ,F12.7,9X, 'L0G(C4)  = '  ,F12 . 7 , 4X , ' (  WEIGHT  =',F8.4,’  )’,//, 

*  T 22 , 'M  =' ,  F  1 2 . 7  ,  /  ) 

RETURN 
END 

SUBROUTINE  PR INT  I  ( N , X , V , G , IER ) 

PRINTING  OF  THE  RESULTS  AND  COMPUTATION/PRINTING  OF  ERRORS 
AND  STANDARD  DEVIATIONS 

DIMENSION  X(9) ,G(9) ,E(900) ,PD(900) ,TC(900) ,W(900) 

COMMON  /P ARM  1 /  NC , ND (5 )  ,  RW( 4 ) 

COMMON  /PARM2/  IC , PLOG ( 900 ) , TLOG ( 900 ) , ULOG ( 900 ) , TAU ( 900 ) , S ( 9 ) 
COMMON  /PARM3/  P ( 900 ) , T ( 900 ) , U ( 900 ) , L ( 20 ) 

COMMON  /PARM4/  P 0 , TO , NDIM , ID ( 5 , 9  ) 

EQUIVALENCE  (E , PLOG ),( PD , TLOG ),( TC , ULOG ) 

EQUIVALENCE  IS  FOR  SPACE  CONSERVATION 

17=0 
TTD=0 . 

TV  =  0 . 

V=SQRT(V /FLOAT (ND(NC+1 ) ) ) 

C 

WRITE (6, 510)  IER, IC,(X(I),G(I), 1=1,5) 

510  FORMAT  (1H0,'»*  RESULTS  OF  COMPUTATION  :  IER  =',I3,4Xt 

*  'SUBROUTINE  FUNCT  WAS  CALLED', 16,'  TIMES  *«', 

*  ///,'  FINAL  VALUES  AND  GRADIENTS 


J 


179 


*  T 35 , ' A  1  =’ ,F12.7,T65, ’ D/D ( A 1 )  =',E15.6,//, 

*T35,'A2  =' ,F12.7,T65, 'D/D(A2)  =',E15.6,//, 

*T35,'A3  =' ,F12.7,T65, ’D/D(A3)  =',E15.6,//, 

*  T35 , ' N  =  ' ,F 12 . 7 , T65 , 'D/D (N )  =',E15.6,//, 

*  T 35 , ' M  =' ,F12.7,T65, ’D/D(M)  =  * , E 1 5 . 6 ) 

WRITE (6,511)  (I,X(I+5) , I,G(I+5) , 1=1 , NC) 

511  FORMAT ( 1  HO, T35, ' LOG ( C • , 1 1 , • )  =  '  , F 1 2 . 7 , T65 , 'D/D(LOG(C' ,11 , 

*  E15.6) 

WRITE ( 6 , 5  1 2  )  V 

512  F0RMAT(/ , 1 H  0 , 'FINAL  STANDARD  DEVIATION  :  ’  ,  F 1 5 . 7  ) 

WRITE ( 6 , 5 1 3 )  C  L ( I ) ,1  =  1 ,20) 

513  FORMAT ( 1 H 0 , T5 , ' (  COMMENT  :  ',20A4,'  )') 

DO  50  J=1 , NC 
J J= J+5 
V  =  0  . 

TD  =  0  . 

IM  =N  D ( J )  +  1 
I N  =  ND ( J  +  1  ) 

K  =  ND ( J  + 1 ) -N D ( J  ) 

RK=FL0AT (K ) 

DO  51  I =IM  ,  IN 

W(I )=X(JJ )+X(4)*PLOG(I )+X(5)*TL0G(I )+UL0G(I) 
R=10.**(X(1)+X(2)*W(I)+X(3)*W(I)**2) 

IFCR.LE.70)  GO  TO  52 
TC ( I )  =  0  . 

GO  TO  53 

52  CONTINUE 

TC ( I  )  =  EXP ( -R  ) 

53  CONTINUE 

E(I )=TAU(I )-TC( I ) 

PD(I)=100.*E(I)/TAU(I ) 

TD=TD+ABS(E(I)) 

V=V+E ( I )  **2 
51  CONTINUE 

TTD=TTD+TD 
TV=TV+V 
TD=TD/RK 
V=V / RK 
SF  =  SQRT ( V ) 

WRITE (6, 514)  (ID(J,I),I=1,9),K,(X(I),I=1,5),J,X(JJ),TD,V,SF 

514  FORMAT  (1H1.T15, 'ACTUAL/COMPUTED  TRANSMITTANCES ',//// , 

*  T  1  0  ,  '  DATA  :  '  ,9A4,5X, '//  OF  POINTS  =  '  ,  14  ,////,  T  1 0  , 

*  ' A  1  ='  , F 1 2 . 7 , 8X ,  'A2  =',F12.7,8X,'A3  =',F12.7,//,T10,'N  =',F12.7 

*  8X  ,  ' M  ='  ,F12.7,8X,  'LOG(C'  ,11  , '  )  =  '  ,  F  1  2 . 7 , //// , 

*  T 1 5 , ’AVERAGE  DISCREPANCY  =',F12.7,//, 

*  T  1  5  ,  'MEAN  SQUARE  ERROR  =',F12.7,//, 

*  T 15 ,  'STANDARD  DEVIATION  ='  ,F12 . 7  ,  ////  , 

*  T  5  ,  '  ,  T 1  1  ,  ’  U  =  '  ,  T26  ,  '  P  =',T36,'T  =',T46,’X  =  '  ,  T56  ,  '  ACTU  A  L  '  , 

*  T65, 'COMPUTED' ,T76, 'DIFFERENCE'  ,T90, '%  DIFF.',/) 

WRITE (6 , 5  1  5  )  (I,U(I),P(I),T(I),W(I),TAU(I),TC(I),E(I),PD(I), 

*  I  =  IM  ,  IN  ) 

515  FORMAT  ( 1 H  , 15 , T10 , El  1 . 4 , T24 , F8 . 2 , T34 , F8 . 2, T43 , F9 . 4 , T53 , F9 . 4 , T63 , 
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*  F9.4,T74,F1 1 . 6 , T87 , F 1 0 . 4 ) 

50  CONTINUE 

FK=FL0AT(ND(NC+1  )  ) 

TTD=TTD/FK 
T V=T  V/FK 
TSD=SQRT (TV) 

WRITE (6 ,516)  NDCNC+1 ) , TTD , TV , TSD 

516  FORMAT  ( 1 H 1 ,///, T 1 0 TOTAL  #  OF  POINTS  USED  = ' , 16 , / / / , T 1 0 , 

*  'GLOBAL  AVERAGE  DISCREPANCY  = ' , F 1 2 . 7 , /// , T 1 0 , 

*  'GLOBAL  MEAN  SQUARE  ERROR  ='  ,  F  1 2 . 7 , /// , T 1 0 , 

*  'GLOBAL  STANDARD  DEVIATION  =',F12.7) 

RETURN 

END 

SUBROUTINE  NMBC ( X , NAME , NC , WN , CF , PO , TO ) 

COMPUTATION  OF  THE  C'-VALUES  FOR  NON-MAJOR  BANDS 

DIMENSION  X(9) ,NAME(20) , WN ( 4  ) 

DIMENSION  CSC  10)  ,FS( 10) 

DF  =  1 . E30 
SGN  =  1  . 

IF ( X ( 3 ) . LT . 0 .  )  SGN  =  - 1  . 

IF  THE  QUADRATIC  TERM  IS  TOO  SMALL,  THEN  IT  WILL  BE  IGNORED 

SMI =-2 . *X ( 3 ) /X (2  ) 

IF(ABSCSMI) . LE . 1 . E-6  )  GO  TO  50 
S  YM  =  1 ./SMI 
50  CONTINUE 

WRITEC6 .5) (NAME ( I ) , 1=1 ,20) 

5  FORMATC1H1 ,T15,20A4) 

WRITE (6 , 10) 

10  FORMAT ( 1 H 0 , T 1 5 ,  '  ***  CALCULATION  OF  THE  SPECTRAL  PARAMETER  FOR', 

*  '  NON-MAJOR  BANDS  *#*',///) 

11  CONTINUE 
NF  REQ  =  0 

12  CONTINUE 

C=0. 

1=0 

15  CONTINUE 

READ (5 ,20, END=40)  KG AS , FREQ , P , T . UGAS , TX 
20  F0RMAT(I2,F10.3,E1 1 . 4 , F9 . 3 , 24X  .  E  1  1  . 4,F7.4) 

Ir(KGAS.EQ.O)  GO  TO  25 
IFCKGAS. LT. 0)  GO  TO  35 
IF(UGAS.GE.DF)  GO  TO  15 


1=1  +  1 
WX=FREQ 
UGAS= JGAS/CF 

IFCSMI . LE . 1 . E-6)  GO  TO  51 

CASE  1  QUADRATIC  TERM  IS  LARGE  AND  USED 
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XS=SYM+SGN«ABS(SQRT(X(2)**2-4.»X(3)«(X(1)-ALOG10(-AL0G(TX)))) 
•  /(2.*X(3) ) ) 

GO  TO  52 

51  CONTINUE 

CASE  2  QUADRATIC  TERM  IS  SMALL  AND  IGNORED 

XS=(ALOG10(-ALOG(TX))-X(1 ))/X(2) 

52  CONTINUE 

XC=X (4 )*AL0G10 (P/P0)+X (5 )*ALOG10 (T0/T)+ALOG10(UGAS) 
C=C+(XS-XC) 

GO  TO  15 

25  C=C/FLOAT(I ) 

NFREQ=NFREQ+1 
CS(NFREQ)=C 
FS(NFREQ)=WX 
DO  27  M=1 , NC 

IF(ABS(WX-WN(M) ) .LE.O. 1 )  CS(NFREQ)=X(5+M) 

27  CONTINUE 

IF(NFREQ.EQ. 10)  GO  TO  30 
GO  TO  12 

30  CONTINUE 

WRITE(6,31)(FS(K) ,K=1 , NFREQ) 

31  FORMAT (1H0.2X, 'WAVE  NUMBER • ,2X , 10F 1 1 . 0) 

WRITE(6,32)(CS(K) ,K=1 , NFREQ) 

32  FORMAT ( 1H0,5X, *C  VALUES  * ,2X, 10F1 1 . 3//) 

GO  TO  11 

35  CONTINUE 

IF(NFREQ.EQ.O)  GO  TO  40 
WRITE(6, 31 ) (FS(K ) ,K=1, NFREQ) 

WRITE (6 , 32) ( CS (K ) ,K=1, NFREQ) 

40  CONTINUE 
RETURN 
END 
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